• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于并行化直接解法的頻率域可控源電磁三維正演

    2015-03-01 01:35:44韓波胡祥云黃一凡彭榮華李建慧蔡建超
    地球物理學(xué)報 2015年8期
    關(guān)鍵詞:場源內(nèi)存電磁

    韓波, 胡祥云, 黃一凡, 彭榮華, 李建慧, 蔡建超

    中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院,地球內(nèi)部多尺度成像湖北省重點實驗室, 武漢 430074

    ?

    基于并行化直接解法的頻率域可控源電磁三維正演

    韓波, 胡祥云*, 黃一凡, 彭榮華, 李建慧, 蔡建超

    中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院,地球內(nèi)部多尺度成像湖北省重點實驗室, 武漢 430074

    電磁法的三維數(shù)值模擬是一個對數(shù)值算法和計算機(jī)硬件要求都非常高的問題.對常用的微分類方法如有限單元法和有限差分法而言,求解最后所得的大型線性方程組是至關(guān)重要的一步,直接影響到正演算法的實用性.如何高效、穩(wěn)定且準(zhǔn)確地解線性方程長期以來一直是被探討的問題.本文實現(xiàn)了基于線性系統(tǒng)直接求解技術(shù)的頻率域可控源電磁(CSEM)三維正演.使用交錯網(wǎng)格有限體積法(FV)來離散化關(guān)于二次電場的Helmholtz方程;使用直接解法取代傳統(tǒng)的迭代解法來求解離散線性系統(tǒng),即對系統(tǒng)矩陣進(jìn)行完全LU分解,具體通過調(diào)用大規(guī)模并行矩陣直接求解器(MUMPS)來實現(xiàn).基于理論模型做了一系列數(shù)值實驗,首先證明了直接解法的高精度和穩(wěn)定性,并考察了其內(nèi)存需求、計算時間和并行可伸縮性等主要計算性能, 最后檢驗了所開發(fā)的算法快速模擬多場源CSEM問題的能力以及對常規(guī)海洋和陸地CSEM模擬的有效性.

    可控源電磁法; 三維模擬; LU分解; 直接解法

    1 引言

    頻率域可控源電磁法(controlled-source electromagnetic methods, CSEM)通過觀測人工場源在地層、空氣和海水等介質(zhì)中激發(fā)的電磁場來獲取介質(zhì)的電性分布,長期以來一直是金屬礦勘探中最重要的地球物理方法(底青云和王若,2008;Hu et al.,2013),并在最近十多年逐漸成為海洋油氣勘探中一種不可或缺的地球物理技術(shù)(Constable,2010).為了通過CSEM觀測數(shù)據(jù)盡可能真實地還原出介質(zhì)的電性結(jié)構(gòu)以提升勘探效果,實施三維勘探并對數(shù)據(jù)進(jìn)行三維模擬是必需的.正演模擬是反演模擬的基礎(chǔ).CSEM中最常用的三維正演方法包括積分方程法(IE)(Avdeev and Knizhnik,2009;王若等,2009;陳桂波等,2009;Zaslavsky et al.,2011),有限單元法(FE)(徐志鋒和吳小平,2010;Schwarzbach et al.,2011;Li et al.,2011;da Silva et al.,2012;Puzyrev et al.,2013),有限差分法(FD)(Newman and Alumbaugh,1995;沈金松,2003;Streich,2009;鄧居智,2011)以及與FD密切相關(guān)的有限體積法(FV)(Weiss and Constable,2006;楊波等,2012;Jahandari and Farquharson,2014;韓波等,2015).對這幾種方法的評述可以在很多文獻(xiàn)中找到(Avdeev,2005;吳桂桔等,2010; B?rner,2010).

    使用FE、FD和FV法對Maxwell方程離散化后一般會得到大型、稀疏線性系統(tǒng)方程,在三維情況下待求解的未知數(shù)很容易達(dá)到上百萬個,要高效地求解這樣的方程并不是一個容易的問題.普遍采取的求解手段是利用Krylov子空間迭代技術(shù)(Saad,2003)來逐步逼近方程的真實解.這一類方法不必顯式存儲系數(shù)矩陣的元素,而只需存儲矩陣與向量的乘積(Weiss,2001),因此對計算機(jī)內(nèi)存的需求很小,完全可以在現(xiàn)代的PC機(jī)上執(zhí)行,并且計算速度很快(在收斂性好的情況下).然而,迭代解法的收斂性能很大程度上取決于系數(shù)矩陣譜的性質(zhì).如果譜的性質(zhì)不佳,就需要大量的迭代來收斂到滿足精度要求的解,但也可能會無法收斂甚至發(fā)散,求解失敗.此時需要對原始方程進(jìn)行"預(yù)優(yōu)"處理來改善系數(shù)矩陣的譜的性質(zhì),再應(yīng)用迭代解法才有可能快速收斂.在電磁法數(shù)值模擬中,空氣、地層和海水之間巨大的電導(dǎo)率差異,以及不同區(qū)域的離散網(wǎng)格尺寸的巨大差異等往往導(dǎo)致最后的線性系統(tǒng)方程譜性質(zhì)極差(比如具有很大的譜條件數(shù)),必須要預(yù)優(yōu).預(yù)優(yōu)因子的選取也是一個問題,好的預(yù)優(yōu)因子往往計算起來較困難,并且沒有一種預(yù)優(yōu)因子能保證迭代解法對于任意復(fù)雜的模型都能穩(wěn)定收斂.

    基于矩陣分解的直接解法是求解線性方程的另一個途徑.與迭代解法相比,直接解法分解矩陣要占用大量的內(nèi)存,且計算時間一般也更長.尤其是對于系統(tǒng)矩陣極其龐大的3D正演,直接解法的這兩個缺點變得非常突出,令普通PC機(jī)難以承受,因此在以往基本不被考慮.然而直接解法相對于迭代解法也有顯著的優(yōu)點,首先是求解精度明顯更高,且計算時間主要取決于系數(shù)矩陣的大小,與其條件數(shù)關(guān)系不大;此外,求解多個系數(shù)矩陣相同而右端項不同的方程組,只需要進(jìn)行一次矩陣分解,因此計算量與解一個方程組相當(dāng)(Grayver,2013),這一點不僅對于多場源的CSEM模擬非常有利(Oldenburg et al.,2008),而且有利于電磁反演中靈敏度的快速構(gòu)建(韓波等,2012).近年來,直接解法在2D/2.5D電磁模擬中已有應(yīng)用(Abubakar et al.,2008;Streich et al.,2011;Luo et al.,2014).隨著數(shù)值算法的發(fā)展,在科學(xué)計算領(lǐng)域出現(xiàn)了許多性能優(yōu)良且免費的矩陣直接求解器軟件包,其中包括針對大規(guī)模稀疏矩陣作了高度優(yōu)化的并行求解器,如MUMPS(Amestoy et al.,2006)和PARDISO(Schenk and G?rtner,2004),加上計算機(jī)性能的快速提升,使得在小型服務(wù)器、工作站和PC機(jī)群上利用直接解法求解3D 電磁問題成為可能.Streich(2009)、da Silva等(2012)以及Jahandari和Farquharson(2014)分別使用FD、FE和FV法進(jìn)行CSEM三維正演計算,他們均利用MUMPS作為線性求解器;Schwarzbach等(2011)則在其高階FE算法中調(diào)用PARDISO解線性方程.

    本文使用交錯網(wǎng)格FV法來進(jìn)行頻率域CSEM的三維正演計算,嘗試?yán)么笠?guī)模并行矩陣求解器MUMPS解線性系統(tǒng)方程.通過大量的理論模型正演測試來考察MUMPS的求解精度、穩(wěn)定性、內(nèi)存需求、計算時間以及并行可伸縮性等計算性能,并驗證所開發(fā)的算法和程序的有效性和廣泛適用性.

    2 方法理論

    2.1 控制方程

    對于一般CSEM勘探所使用的頻率范圍和大部分已知的傳播介質(zhì)的電阻率變化范圍,位移電流的影響可以忽略,則在電流型場源激發(fā)下關(guān)于電場的二階頻率域微分方程為(取時諧因子為eiω t)

    (1)

    其中,E為電場強(qiáng)度(V/m),σ為電導(dǎo)率(S/m),ω為角頻率,μ0為真空中的磁導(dǎo)率(H/m),J為外加場源的電流密度(A·m-2).J為沖激函數(shù),在場源點具有奇異性(無窮大)(WeissandConstable,2006),使得數(shù)值求解方程(1)有一定的困難.為此,可假設(shè)有一個參考模型σp,在場源J的激發(fā)下同樣滿足上述形式的方程:

    (2)

    式(1)與式(2)相減,可得

    (3)

    若將實際模型的場看做總場,參考模型的場看做一次場,令二次場Es=E-Ep,式(3)可重寫為

    (4)

    其中Js=(σ-σp)Ep,可看做激發(fā)二次場的電流密度.方程(4)便是關(guān)于二次電場的二階控制方程,也是本文使用FV法直接求解的方程.與總場方程(1)相比,由于Js不是奇異函數(shù),方程(4)求解起來更加容易.對于一次場,一般將參考模型選為所需計算模型的一維“背景模型”,可以通過解析法求得;一次場與二次場相加便得到總場.

    2.2 有限體積法

    圖1 (a)電磁場的交錯采樣和(b)采樣點處的控制體積Fig.1 (a) Sampling positions of EM fields on a staggered grid and (b) the control volume for a sampling site

    求解形如式(4)的Helmholtz方程,一種簡單而有效的方法是交錯網(wǎng)格有限差分法(Staggered Finite-Difference, SFD),其將研究區(qū)域離散化為一系列規(guī)則長方體單元,然后在長方體單元上交錯地對電場和磁場采樣,如圖1a所示.該方法在電磁法3D正演中被廣泛使用,其基本理論比較簡單,可以在很多文獻(xiàn)中找到(e.g.,Newman and Alumbaugh,1995),這里不再贅述.使用SFD法對(4)進(jìn)行離散化直接得到的線性方程系數(shù)矩陣是不對稱的,為了使其對稱以便于方程的求解,每個采樣點的線性方程兩端需要同時乘以一個對應(yīng)于該采樣點的“尺度因子”V,V具有某種體積的意義,即SFD線性方程實際上對應(yīng)的微分方程為

    (5)

    對SFD法稍加改動,便可以實現(xiàn)交錯網(wǎng)格有限體積法(SFV),即在每個采樣點處定義一個包圍該采樣點的控制體積,由共享該采樣點所在棱邊的相鄰4個網(wǎng)格單元體積的1/4組成.以x分量為例,如圖1b所示.對式(4)在某個體積單元Ω內(nèi)求積分:

    (6)

    然后對式(6)求解.事實上,式(5)中的尺度因子剛好是SFV中的控制體積值,因此式(5)與式(6)在物理意義上的區(qū)別在于前者用某個離散點處的變量值代表包圍該點的一片區(qū)域內(nèi)每一處的變量值,后者則考慮到變量在區(qū)域內(nèi)的變化,對區(qū)域進(jìn)行體積分.式(6)左端第一項可以完全展開(楊波等,2012),并且與(5)左端第一項具有完全相同的離散形式.對于(6)的左端第二項,通常認(rèn)為二次電場在控制體積范圍內(nèi)的變化不大,因此Es可提取到積分號以外;采樣點處的電導(dǎo)率可與SFD中一樣,使用采樣點所在棱邊的相鄰4個單元電導(dǎo)率的體積加權(quán)平均值,也不必求積分.經(jīng)過這樣的簡化,SFV的左端項就與SFD完全相同,因此二者的線性方程系數(shù)矩陣相同,唯一不同的是右端項.由于CSEM的電磁場可能在很小的空間范圍內(nèi)發(fā)生較大變化(如場源附近),因此SFV比SFD更能準(zhǔn)確地代表電磁場的分布.

    CSEM的場源是局部的,因此可讓剖分域的邊界離場源足夠遠(yuǎn),從而運用齊次的Dirichlet邊界條件,即令邊界上的切向電場為0;SFV法右端項的一次場體積積分可使用數(shù)值積分實現(xiàn),例如Gauss-Legendre積分(DavisandRabinowitz,2007);一次場的計算為一維正演問題,使用基于Hankel變換的準(zhǔn)解析法(Key,2009).

    2.3 大型線性系統(tǒng)方程的求解

    Helmholtz方程經(jīng)SFD或SFV離散化后得到線性系統(tǒng)方程:

    Ae=b,

    (7)

    其中A為大型、對稱、稀疏復(fù)矩陣(僅對角元為復(fù)數(shù)),其行(列)數(shù)也即待求的未知數(shù)個數(shù)約為3Nx·Ny·Nz(Nx、Ny和Nz分別為三個方向的網(wǎng)格數(shù)),每行最多只有13個非零元.作為一個簡單示例,圖2中顯示了一個網(wǎng)格數(shù)為5×6×7的模型的SFD/SFV離散線性方程系數(shù)矩陣的稀疏結(jié)構(gòu).

    本文使用基于矩陣三角分解的直接解法求解方程(7),具體通過一個開源的軟件包——多波前大規(guī)模并行求解器(MultifrontalMassivelyParallelSolver,MUMPS)(Amestoyetal.,2006)來實現(xiàn).MUMPS專門針對大型稀疏矩陣,并且還能利用矩陣的對稱性,就內(nèi)存需求和運行時間作了高度優(yōu)化,因此非常適合解方程(7);此外MUMPS還能通過消息傳遞接口(MPI)實現(xiàn)并行化求解.MUMPS解方程包括3個階段:①矩陣的分析與預(yù)處理;②矩陣的LU(LDLT)分解;③三角方程的求解.其中矩陣的分解占用絕大部分計算時間和內(nèi)存.

    圖2 網(wǎng)格數(shù)為5×6×7的模型的SFD/SFV離散線性方程系數(shù)矩陣A的稀疏結(jié)構(gòu)Fig.2 Sparsity pattern of A arising from the SFD/SFV discretization on a staggered-grid with size 5×6×7

    對照控制方程(4)或(6)來看方程(7),在模型剖分不變的情況下,改變頻率,則系數(shù)矩陣A與右端項b都會改變,需要重新分解A并求解;若模型剖分和頻率均固定,只改變發(fā)射源的方位或位置,則一次場會改變,導(dǎo)致b發(fā)生改變,但A不會受到影響,因此無需重新進(jìn)行矩陣分解,而可以利用原有的分解因子,對新的右端項重新進(jìn)行解三角方程就夠了.這樣就大大節(jié)省了計算量.在實際CSEM勘探中(Oldenburgetal.,2008;Constable,2010;Grayver,2013),為了更全面地獲取地下信息,常需要頻繁改變發(fā)射源的極化方向(方位)和空間位置,此時利用直接解法可以快速計算多方位、多空間位置發(fā)射源的電磁響應(yīng).

    2.4 計算流程

    根據(jù)前文的描述,本文的SFV算法程序的計算流程如圖3所示,從中可以清晰地看出MUMPS的各個階段需要執(zhí)行的次數(shù)與頻率數(shù)、場源數(shù)之間的關(guān)系.

    圖3 本文SFV三維正演算法的計算流程Fig.3 Workflow of the presented SFV modeling scheme

    3 直接解法的性能測試

    為了測試MUMPS求解大型線性系統(tǒng)的性能,基于兩個較典型的層狀地電模型(圖4)進(jìn)行了一系列正演測試.圖4a為海洋CSEM中常用的模擬高阻油氣薄層的模型(Constable and Weiss,2006),圖4b為一個陸地三層H型模型.本節(jié)的測試主要與模型參數(shù)和頻率有關(guān),而與一些具體的收發(fā)設(shè)置如測點設(shè)置無關(guān),因此將不提及這些設(shè)置.本文的所有計算均在中科曙光“天闊W580I-G10”工作站上完成,該工作站包含兩顆Intel XEON E5-2603型CPU,每顆CPU為4核心4線程,主頻1.8 GHz.為了應(yīng)對矩陣分解的巨大內(nèi)存需求,配置了128 GB的內(nèi)存.操作系統(tǒng)為Ubuntu 12.04.

    3.1 解方程的精度與穩(wěn)定性

    從基本原理上看,直接解法未作任何近似,因此其求解結(jié)果在理論上是精確的,但實際執(zhí)行時由于計算機(jī)的字長有限,對浮點數(shù)運算會有舍入誤差,導(dǎo)致直接解法也不可避免地存在誤差,但這種誤差與迭代解法中由于對數(shù)學(xué)問題本身的近似帶來的誤差有本質(zhì)區(qū)別.

    圖4 用來測試直接解法性能的(a)海洋與(b)陸地一維層狀模型Fig.4 The (a) marine layered 1D model and (b) land layered 1D model used to test the performance of the direct solver

    方程(7)中,系數(shù)矩陣A的數(shù)值特性與頻率、空氣電導(dǎo)率等密切相關(guān),因此方程的求解也可能會受這些因素的影響,尤其是迭代解法往往對這些因素非常敏感.為此,對以上兩個1D模型正演時,大范圍改變頻率和空氣電導(dǎo)率來測試直接解法的穩(wěn)定性.令空氣電導(dǎo)率從10-5S/m變化到10-11S/m,頻率在實際勘探常用的范圍內(nèi)變化,表1和表2給出了直接解法解方程的殘差和矩陣分解耗時.這里殘差的定義式為‖Ae-b‖,相對殘差為‖Ae-b‖/‖b‖;兩個模型的網(wǎng)格剖分?jǐn)?shù)分別為68×40×43和64×42×46;使用了8個線程并行計算.可以看出,空氣電導(dǎo)率的變化對相對殘差的影響非常?。活l率由高到低變化,相對殘差會有明顯的增加,但最大的相對殘差也僅僅在10-11數(shù)量級,而迭代解法收斂的相對殘差門檻值一般設(shè)在10-5~10-8之間(Grayver,2013).在矩陣分解耗時方面,頻率和空氣電導(dǎo)率變化對海洋1D模型的影響很小,對陸地模型的影響稍大,但時間的波動相對于整體耗時并不大.由此可見,在實際勘探常用的頻率范圍內(nèi),對于包含很大電性差異的模型,直接解法都具有相當(dāng)高的精度,計算速度也很穩(wěn)定.

    表1 對海洋1D模型,頻率與空氣層電導(dǎo)率在較大范圍內(nèi)變化時直接解法的精度和矩陣分解時間Table 1 The accuracy and factorization time of the direct solver for the marine 1D model for a wide range of frequencies & air conductivities

    3.2 內(nèi)存需求與計算時間

    再來考察直接解法的內(nèi)存需求與計算時間對模型剖分網(wǎng)格數(shù)的敏感程度.對海洋1D模型,發(fā)射頻率為1 Hz,空氣電導(dǎo)率取10-11S/m(本文從此處開始,空氣電導(dǎo)率值均取該值).設(shè)計了從小到大21種網(wǎng)格,其中最小的網(wǎng)格數(shù)為7×8×8,最大的網(wǎng)格數(shù)為98×80×57.圖5顯示了直接解法的內(nèi)存需求與計算時間(MUMPS的三個階段之和)隨未知數(shù)個數(shù)(約為3Nx·Ny·Nz)的變化.在調(diào)用MUMPS時執(zhí)行了“Out-Of-Core”選項,即利用硬盤來存儲矩陣的分解因子,這樣能節(jié)約很大一部分內(nèi)存.對于較大的網(wǎng)格數(shù),使用了多個進(jìn)程并行計算(進(jìn)程的增加會導(dǎo)致內(nèi)存需求的增加,見下一小節(jié)).可以看到,內(nèi)存需求與計算時間均隨模型網(wǎng)格數(shù)的增加而快速增加.當(dāng)網(wǎng)格數(shù)為78×50×57時(對應(yīng)的未知數(shù)個數(shù)約為650000),內(nèi)存需求在16GB左右,而當(dāng)網(wǎng)格數(shù)達(dá)到98×80×57時,內(nèi)存需求則達(dá)到了接近120GB.雖然無法給出網(wǎng)格數(shù)與內(nèi)存需求之間的具體關(guān)系式,但不難看出整體上內(nèi)存需求增加的速度要遠(yuǎn)快于網(wǎng)格數(shù)增加的速度.

    表2 對陸地1D模型,頻率與空氣層電導(dǎo)率在較大范圍內(nèi)變化時直接解法的精度和矩陣分解時間Table 2 The accuracy and factorization time of the direct solver for the land 1D model for a wide range of frequencies & air conductivities

    圖5 對海洋1D模型,對于不同的模型網(wǎng)格剖分?jǐn)?shù)直接解法的(a)內(nèi)存需求與(b)解方程時間Fig.5 (a) Memory requirement and (b) solving time of the direct solver for the marine 1D model for various grid sizes

    3.3 并行可伸縮性

    在并行計算中,可伸縮性(Scalability)是指系統(tǒng)通過改變可用計算資源和調(diào)度方式來動態(tài)調(diào)整自身計算性能的能力.為了測試MUMPS的可伸縮性,依次使用1~8個進(jìn)程對海洋1D模型進(jìn)行正演計算,頻率為1 Hz,使用3種網(wǎng)格剖分,統(tǒng)計了解方程的時間、相對于單進(jìn)程運算的加速因子、總內(nèi)存消耗和平均每個進(jìn)程的內(nèi)存消耗等,如圖6所示.可以看到,計算時間隨進(jìn)程數(shù)的增加穩(wěn)定下降,當(dāng)進(jìn)程數(shù)達(dá)到8時,計算時間是單進(jìn)程的1/4左右.但從加速因子來看,還遠(yuǎn)不能達(dá)到理論上最理想的并行(加速因子等于進(jìn)程數(shù)),這是由矩陣分解的可并行程度有限所決定的.對于同一種網(wǎng)格,總體內(nèi)存消耗隨進(jìn)程數(shù)的增加呈上升趨勢(個別情況除外),8進(jìn)程運行的總內(nèi)存消耗是單進(jìn)程的4倍左右,這是由于進(jìn)程之間需要共享很多變量.因此在試圖增加進(jìn)程數(shù)來提升計算速度時,還要考慮到進(jìn)程增多帶來的內(nèi)存需求增大.盡管如此,內(nèi)存消耗平均到每個進(jìn)程上,還是會隨進(jìn)程數(shù)增加而下降的,這一點對機(jī)群系統(tǒng)比較有利.

    4 三維正演程序的有效性檢驗

    接下來通過具體正演結(jié)果的比較與分析來檢驗所開發(fā)的SFV三維正演程序的有效性.

    4.1 與解析解的對比

    首先,還是基于第3節(jié)中的兩個層狀模型,比較SFV數(shù)值解與解析解.所有的解析解均使用開源的Dipole1D程序(Key,2009)獲得.

    對于海洋層狀模型,觀測系統(tǒng)設(shè)置為:發(fā)射電偶極沿x方向,位于海底上方50 m,取偶極中心在x-y平面內(nèi)的投影為坐標(biāo)原點,發(fā)射頻率為1 Hz;接收器位于海底表面,沿x軸方向布設(shè)(即inline觀測,參考圖4a).圖7給出了網(wǎng)格數(shù)為86×50×57時Ex分量和Hy分量的振幅-收發(fā)距(MVO)曲線與相位-收發(fā)距(PVO)曲線的SFV解與解析解對比,圖8則為在3種不同網(wǎng)格剖分下SFV解相對于解析解的誤差.可以看出,SFV解的誤差整體上隨收發(fā)距的增加而增大;網(wǎng)格為54×32×27時誤差非常明顯,Hy分量振幅的最大相對誤差可達(dá)20%以上;當(dāng)網(wǎng)格加密為68×40×43時,數(shù)值解的精度得到了極大的提高,振幅最大相對誤差在8%左右,相位最大相差6°左右;進(jìn)一步細(xì)化網(wǎng)格至86×50×57時,精度繼續(xù)得以提升,但并不顯著,這也許跟網(wǎng)格加密的方式有一定關(guān)系(不在本文討論范圍以內(nèi)).

    圖6 對海洋層狀模型,直接解法解方程的時間與內(nèi)存需求隨所使用的進(jìn)程數(shù)的改變Fig.6 Solving time & memory requirement of the direct solver vary as the number of computing processors changes for the marine layered model

    圖7 對海洋層狀模型,網(wǎng)格為86×50×57時水平(a, b)電場和(c, d)磁場的SFV數(shù)值解(圓圈)與解析解(線條)對比Fig.7 Comparison of horizontal (a, b) electric fields and (c, d) magnetic fields obtained from the SFV method (circles) and the analytic solution (line) for the marine 1D model with grid size 86×50×57

    圖8 對于海洋層狀模型,在3種不同網(wǎng)格剖分下SFV解相對于解析解的誤差Fig.8 Relative errors of SFV solutions compared to analytic solutions for the marine 1D model for three different grid sizes

    圖9 對陸地層狀模型,坐標(biāo)為(8000, 0)的測點的水平(a, b)電場、(c, d)磁場和(e, f)卡尼亞視電阻率的SFV解(圓圈)與解析解(線條)對比,以及(g, h)視電阻率的SFV解相對于解析解的誤差Fig.9 Comparison of (a, b) electric fields, (c, d) magnetic fields and (e, f) Cagniard apparent resistivities of the site with coordinate (8000, 0) obtained from the SFV method (circles) and the analytic solution (line) for the land 1D model. (g, h) are the apparent resistivity errors of SFV solutions compared to analytic solutions

    對于陸地層狀模型,采用地面可控源音頻大地電磁法(CSAMT)(底青云和王若,2008)的觀測方式,觀測系統(tǒng)設(shè)置為:發(fā)射源為沿y方向的1000 m的接地導(dǎo)線,取其中心為坐標(biāo)原點,發(fā)射0.125~8192 Hz之間按對數(shù)等間隔分布的17個頻率;測線平行于發(fā)射源布設(shè),中點到場源中心的最大距離為8 km(參考圖4b).網(wǎng)格數(shù)為66×80×46.圖9給出了坐標(biāo)為(8000, 0)的測點上所有頻率的Ey分量、Hx分量和卡尼亞視電阻率ρyx的SFV解與解析解對比,以及ρyx的誤差.其中ρyx的計算式為

    (8)

    可以看出,數(shù)值解與解析解吻合較好.頻率位于64~256 Hz之間時誤差稍大,其余頻率的視電阻率相對誤差均在8%以內(nèi),阻抗相位誤差均在±2°以內(nèi).

    4.2 多場源問題的模擬

    考慮海洋CSEM勘探中需要發(fā)射源在多個位置激發(fā)的情況,對上述海洋1D模型,令發(fā)射源沿x方向在-2000 m到8000 m之間移動,其余觀測參數(shù)與4.1節(jié)中相同.如2.3節(jié)中描述的,頻率不變時,發(fā)射源每改變一次激發(fā)的位置(或方位),就會給三維正演模擬中線性系統(tǒng)方程增加一個右端項,但不會增加系數(shù)矩陣.圖10為右端項數(shù)目由1增加到101時,SFV正演計算的總時間及各主要部分消耗的時間變化.可以看到,當(dāng)右端項比較少時,矩陣分解占用了絕大部分計算時間,而解三角方程(Solution)這一步消耗的時間幾乎可以忽略不計;隨著右端項的增多,由于SFV法每構(gòu)建一次右端項都要在控制體積內(nèi)進(jìn)行數(shù)值積分,構(gòu)建右端項所占用的時間比例越來越大;當(dāng)右端項增加到101個時,構(gòu)建右端項占用的時間比例已經(jīng)很顯著了,但解三角方程所占比重仍然很小,總計算時間也僅僅只比1個右端項時增加了約1倍.

    4.3 海洋三維模型

    若將前面的海洋層狀模型中高阻層在水平方向的邊界截斷,則得到更接近真實情況的三維油氣儲層模型,如圖11所示.已有不少其他作者研究過類似的模型(e.g.,Constable and Weiss,2006;Streich,2009),本文的模型與這些作者略有差別,即高阻體在x-y平面內(nèi)的截面為正方形而非圓形.

    分別考慮海水無窮深、不含空氣層和海水1000 m深、包含空氣層這兩種情況.參數(shù)設(shè)置為:場源沿x方向,位于高阻體x負(fù)向邊界中點的正上方,海底上方50 m,坐標(biāo)為(0, 0, -50 m);發(fā)射頻率為1 Hz;接收器沿x方向布設(shè)、分布于海底y=0的測線上(inline觀測).對高阻體截面邊長(設(shè)為a)為1 km、2 km和5 km的三種情況進(jìn)行了正演.所得到的水平電場Ex分量如圖12所示,其中還給出了不含高阻體的半空間(a=0)和高阻體在水平方向無限延伸(a=∞)這兩種極限情況的響應(yīng)(使用解析法獲得).

    看圖12a和12b,Ex響應(yīng)的分段特征非常明顯:在小收發(fā)距內(nèi),經(jīng)海水傳播的直達(dá)波占主導(dǎo)地位,振幅大約按1/r3衰減;隨著收發(fā)距的增大,地層波開始起主導(dǎo)作用,含高阻層(a>0)與不含高阻層(a=0)的響應(yīng)逐漸區(qū)分開來,振幅曲線在對數(shù)坐標(biāo)系中的斜率趨于常數(shù),說明振幅趨于按指數(shù)衰減.a=1 km與半空間模型的響應(yīng)幾乎重合,說明該尺寸的異常體難以分辨;a=2 km與半空間的響應(yīng)雖然在大收發(fā)距時的衰減速度保持一致,但二者的絕對差值較為明顯,說明該尺寸的異常體容易分辨;a=5 km與半空間的響應(yīng)差別非常明顯,且在收發(fā)距大于高阻體邊長時與高阻體無限延伸(a=∞)的響應(yīng)區(qū)分開來,衰減速度逐漸趨于半空間模型的衰減速度.根據(jù)以上分析,可以再次確認(rèn)Constable和Weiss

    圖10 對于海洋層狀模型,SFV正演總體計算時間及各主要階段消耗的時間隨右端項數(shù)目的變化Fig.10 The total running time and running time of the significant tasks in SFV modeling vs the numbers of right hand sides for the marine 1D model

    (2006)的結(jié)論:對于此類水平薄板狀高阻體,只有當(dāng)其水平尺寸至少為埋深的2倍時,才可以利用CSEM將其分辨出來.

    再來看包含空氣層的情況.對比圖12c、d和a、b,可以看到,空氣波的影響非常明顯,在收發(fā)距超過5 km時,a=1 km、a=2 km和半空間模型的響應(yīng)中空氣波開始起主導(dǎo)作用,電場振幅衰減極其緩慢;而a=5 km的模型直到收發(fā)距達(dá)到8 km左右時才開受空氣波影響,a=∞的模型在所展示的收發(fā)距內(nèi)還未被空氣波影響.

    圖11 海洋三維油氣儲層模型Fig.11 The marine 3D reservoir model

    圖12 三維儲層模型的儲層截面邊長變化時的Ex響應(yīng)(a, b)為不考慮空氣層、假設(shè)海水無限深的情況,(c, d) 為海水1000 m深、包含空氣層的情況.Fig.12 Ex responses for the 3D reservoir model with various horizontal edge lengths of the reservoir(a, b) is for the case air layer is not presented and the seawater is arbitrarily thick, while (c, d) is for the case the seawater is 1000 m thick and air layer is included.

    4.4 陸地三維模型

    最后對一個陸地三維模型進(jìn)行正演.如圖13所示,在100 Ωm的均勻半空間中植入一個400 m×500 m×400 m的低阻塊體,埋深為200 m,電阻率5 Ωm.使用兩個分別沿y方向和x方向的正交場源交替發(fā)射來模擬張量CSAMT觀測,發(fā)射導(dǎo)線長度均為1000 m,發(fā)射0.125~8192 Hz之間按對數(shù)等間隔分布的17個頻率;發(fā)射源的中心為坐標(biāo)原點,到低阻體中心的水平距離為6000 m;網(wǎng)格數(shù)為62×44×53.

    由于有兩個正交場源,因此可以使用類似于大地電磁法(MT)中的方式來計算張量阻抗,進(jìn)而計算卡尼亞視電阻率響應(yīng),具體計算式見譚捍東等(2004).圖14和圖15分別為沿x方向和沿y方向穿過低阻體正上方的兩條測線上的響應(yīng),圖16則為頻率為4 Hz的響應(yīng)切片,均顯示了XY和YX兩種模式的響應(yīng).整體上來看,這些對CSEM的電磁場用MT的方式處理得到的響應(yīng)與MT響應(yīng)有一定的相似度,異常體的空間輪廓在正演響應(yīng)中有良好的反映.XY和YX模式的響應(yīng)特征有所不同,XY響應(yīng)對x方向的電性邊界反映較好,YX響應(yīng)則對y方向的電性邊界反映較好;在與低阻體對應(yīng)的位置、對應(yīng)的頻率處,XY視電阻率比YX視電阻率在數(shù)值上更接近低阻體的電阻率值.在頻率低于2 Hz時,視電阻率值均超過了300 Ωm,甚至可達(dá)1000 Ωm,這是因為頻率較低時,電磁波在背景地層中的趨膚深度較大,測點所處的位置遠(yuǎn)遠(yuǎn)不能滿足遠(yuǎn)區(qū)條件,此時卡尼亞視電阻率公式完全不適用.

    圖13 陸地三維低阻體模型以及發(fā)射源的布置示意圖Fig.13 The land 3D conductive brick model and the source geometry

    5 結(jié)論與討論

    本文利用多波前大規(guī)模并行求解器軟件包MUMPS,實現(xiàn)了基于矩陣三角分解的頻率域可控源電磁三維有限體積正演.

    針對MUMPS直接解法的計算性能做了一系列測試,結(jié)果表明:(1)直接解法非常穩(wěn)定,不易被系數(shù)矩陣的病態(tài)程度所影響,在實際CSEM勘探常用的頻率范圍內(nèi),對于包含很大電性差異的模型,都能獲得高精度的解,計算速度也很穩(wěn)定;(2)直接解法的內(nèi)存需求與計算時間取決于模型網(wǎng)格數(shù),且均隨模型網(wǎng)格數(shù)的增加而急劇增加,內(nèi)存需求很容易超出單臺普通PC機(jī)的物理內(nèi)存;(3)使用多進(jìn)程并行執(zhí)行MUMPS時,其計算速度能得到大幅度提升,也能有效地減少單節(jié)點的內(nèi)存消耗,但由于矩陣分解本身的可并行程度有限,MUMPS的并行可伸縮性也比較有限.

    對理論模型正演結(jié)果的比較與分析表明:本文的SFV三維正演算法是有效的,適用于常規(guī)的陸地和海洋CSEM勘探,并且由于利用了直接解法,非常適合于CSEM中典型的多場源問題,場源的增加只會帶來極小的額外計算量.

    就目前最普通、最可用的計算設(shè)備如個人電腦、工作站和小型服務(wù)器的平均性能而言,電磁法的三維數(shù)值模擬仍然是一個對計算資源需求很高的問題,使用直接解法時更是如此,因此模型剖分的網(wǎng)格數(shù)要嚴(yán)格控制(本文3.2節(jié)中對此給出了參考),這無疑會影響計算結(jié)果的準(zhǔn)確性.本文的研究便受限于計算設(shè)備的性能,在并行計算的進(jìn)程數(shù)、地面CSAMT模擬中的收發(fā)距等方面受到了較大的限制.另外,在計算設(shè)備的性能足夠強(qiáng)勁的情況下,對多頻率的正演計算實施“多層次”并行化無疑能進(jìn)一步提高效率,即將不同頻率的正演任務(wù)分配到不同的節(jié)點組,然后每個節(jié)點組對單個頻率進(jìn)行并行正演計算.這些都是有待改進(jìn)和完善的地方.

    隨著數(shù)值算法的發(fā)展以及計算機(jī)性能的快速提升,可以預(yù)見,直接解法由于其特殊的優(yōu)點,必將越來越多地被應(yīng)用于3D電磁模擬中.

    致謝 感謝MUMPS軟件包和Dipole1D的開發(fā)者們,沒有他們的開源精神,本文的研究將無法完成;特別感謝美國俄勒岡州立大學(xué)(OSU)Adam Schultz教授的指導(dǎo).

    圖14 陸地三維模型在x=6000 m的測線上的卡尼亞視電阻率響應(yīng)擬斷面圖Fig.14 The Cagniard apparent resistivity response pseudo-sections of the survey line at x=6000 m for the land 3D model

    圖15 陸地三維模型在y=0 m的測線上的卡尼亞視電阻率響應(yīng)擬斷面圖Fig.15 The Cagniard apparent resistivity response pseudo-sections of the survey line at y=0 m for the land 3D model

    圖16 陸地三維模型在頻率f=4 Hz時的卡尼亞視電阻率響應(yīng)頻率切片圖Fig.16 The Cagniard apparent resistivity response frequency slices at f=4 Hz for the land 3D model

    Abubakar A, Habashy T M, Druskin V L, et al. 2008. 2.5D forward and inverse modeling for interpreting low-frequency electromagnetic measurements.Geophysics, 73(4): F165-F177.

    Amestoy P R, Guermouche A, L′Excellent J, et al. 2006. Hybrid scheduling for the parallel solution of linear systems.ParallelComputing, 32(2): 136-156.

    Avdeev D B. 2005. Three-dimensional electromagnetic modelling and inversion from theory to application.SurveysinGeophysics, 26(6): 767-799. Avdeev D, Knizhnik S. 2009. 3D integral equation modeling with a linear dependence on dimensions.Geophysics, 74(5): F89-F94. B?rner R U. 2010. Numerical modelling in Geo-Electromagnetics: advances and challenges.SurveysinGeophysics, 31(2): 225-245. Chen G B, Wang H N, Yao J J, et al. 2009. Modeling of electromagnetic responses of 3-D electrical anomalous body in a layered anisotropic earth using integral equations.ChineseJ.Geophys. (in Chinese), 52(8): 2174-2181, doi: 10.3969/j.issn.0001-5733.2009.08.028.

    Constable S C, Weiss C J. 2006. Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling.Geophysics, 71(2): G43-G51.

    Constable S C. 2010. Ten years of marine CSEM for hydrocarbon exploration.Geophysics, 75(5): 75A67-75A81.

    da Silva N V, Morgan J V, MacGregor L, et al. 2012. A finite element multifrontal method for 3D CSEM modeling in the frequency domain.Geophysics, 77(2): E101-E115.

    Davis P J, Rabinowitz P. 2007. Methods of Numerical Integration. 2nd ed. New York: Dover Publications.

    Deng J Z. 2011. Studies of three dimensional stagged-grid finite difference for controlled source audio-frequency magnetotelluric numerical simulation [Ph. D. thesis] (in Chinese). Beijing: China University of Geosciences (Beijing).

    Di Q Y, Wang R. 2008. CSAMT Forward Modeling and Inversion and Its Application (in Chinese). Beijing: Science Press. Grayver A V. 2013. Three-dimensional controlled-source electromagnetic inversion using modern computational concepts [Ph. D. thesis]. Berlin: Free University of Berlin.

    Han B, Hu X Y, He Z X, et al. 2012. Mathematical classification of magnetotelluric inversion methods.OilGeophysicalProspecting(in Chinese), 47(1): 177-187.

    Han B, Hu X Y, Schultz A, et al. 2015. Three-dimensional forward modeling of the marine controlled-source electromagnetic field with complex source geometries.ChineseJ.Geophys. (in Chinese), 58(3): 1059-1071, doi: 10.6038/cjg20150330.

    Hu X Y, Peng R H, Wu G J, et al. 2013. Mineral Exploration using CSAMT data: Application to Longmen region metallogenic belt, Guangdong Province, China.Geophysics, 78(3): B111-B119. Jahandari H, Farquharson C G. 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids.Geophysics, 79(6): E287-E302. Key K. 2009. 1D inversion of multicomponent, multifrequency marine CSEM data: methodology and synthetic studies for resolving thin resistive layers.Geophysics, 74(2): F9-F20.

    Li J H, Zhu Z Q, Liu S C, et al. 2011. 3D numerical simulation for the transient electromagnetic field excited by the central loop based on the vector finite-element method.JournalofGeophysicsandEngineering, 8(4): 560-567. Luo M, Li Y, Liu Y. 2014. 2.5D controlled-source electromagnetic modeling by an adaptive finite-element method with an arbitrarily oriented finite length dipole source. 22nd EM Induction Workshop. Abstracts. Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences.GeophysicalProspecting, 43(8): 1021-1042.

    Oldenburg D W, Haber E, Shekhtman R. 2008. Forward modelling and inversion of multi-source TEM data. ∥ 78th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 559-563. Puzyrev V, Koldan J, de la Puente J, et al. 2013. A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling.Geophys.J.Int., 193(2): 678-693. Saad Y. 2003. Iterative Methods for Sparse Linear Systems. 2nd ed. Philadelphia: SIAM.

    Schenk O, G?rtner K. 2004. Solving unsymmetric sparse systems of linear equations with PARDISO.FutureGenerationComputerSystems, 20(3): 475-487.

    Schwarzbach C, B?rner R U, Spitzer K. 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—a marine CSEM example.Geophys.J.Int., 187(1): 63-74. Shen J S. 2003. Modeling of 3-D electromagnetic responses in frequency domain by using staggered-grid finite difference method.ChineseJ.Geophys. (in Chinese), 46(2): 281-289.

    Streich R. 2009. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy.Geophysics, 74(5): F95-F105.

    Streich R, Becken M, Ritter O. 2011. 2.5D controlled-source EM modeling with general 3D source geometries.Geophysics, 76(6): F387-F393.

    Tan H D, Wei W B, Deng M, et al. 2004. General use formula in MT tensor impedance.OilGeophysicalProspecting(in Chinese), 39(1): 113-116. Wang R, Di Q Y, Wang M Y, et al. 2009. Research on the effect of 3D body between transmitter and receivers on CSAMT response using Integral Equation method.ChineseJ.Geophys. (in Chinese), 52(6): 1573-15821, doi: 10.3969/j.issn.0001-5733.2009.06.019.

    Weiss C J. 2001. A matrix-free approach to solving the fully 3D electromagnetic induction problem. ∥71th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1451-1454. Weiss C J, Constable S C. 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part II — Modeling and analysis in 3D.Geophysics, 71(6): G321-G332.

    Wu G J, Hu X Y, Liu H. 2010. Progress in CSAMT three-dimensional forward numerical simulation.ProgressinGeophys. (in Chinese), 25(5): 1795-1801, doi: 10.3969/j.issn.1004-2903.2010.05.037.

    Xu Z F, Wu X P. 2010. Controlled source electromagnetic 3-D modeling in frequency domain by finite element method.ChineseJ.Geophys. (in Chinese), 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019.

    Yang B, Xu Y X, He Z X, et al. 2012. 3D frequency-domain modeling of marine controlled source electromagnetic responses with topography using finite volume method.ChineseJ.Geophys. (in Chinese), 55(4): 1390-1399, doi: 10.6038/j.issn.0001-5733.2012.04.035.

    Zaslavsky M, Druskin V, Davydycheva S, et al. 2011. Hybrid finite-difference integral equation solver for 3D frequency domain anisotropic electromagnetic problems.Geophysics, 76(2): F123-F137.

    附中文參考文獻(xiàn)

    陳桂波, 汪宏年, 姚敬金等. 2009. 用積分方程法模擬各向異性地層中三維電性異常體的電磁響應(yīng). 地球物理學(xué)報, 52(8): 2174-2181, doi: 10.3969/j.issn.0001-5733.2009.08.028.

    鄧居智. 2011. 可控源音頻大地電磁法三維交錯采樣有限差分?jǐn)?shù)值模擬研究[博士論文]. 北京: 中國地質(zhì)大學(xué)(北京).

    底青云, 王若. 2008. 可控源音頻大地電磁數(shù)據(jù)正反演及方法應(yīng)用. 北京: 科學(xué)出版社.

    韓波, 胡祥云, 何展翔等. 2012. 大地電磁反演方法的數(shù)學(xué)分類. 石油地球物理勘探, 47(1): 177-187.

    韓波, 胡祥云, Schultz A等. 2015. 復(fù)雜場源形態(tài)的海洋可控源電磁三維正演. 地球物理學(xué)報, 58(3): 1059-1071, doi: 10.6038/cjg20150330.

    沈金松. 2003. 用交錯網(wǎng)格有限差分法計算三維頻率域電磁響應(yīng). 地球物理學(xué)報, 46(2): 281-289.

    譚捍東, 魏文博, 鄧明等. 2004. 大地電磁法張量阻抗通用計算公式. 石油地球物理勘探, 39(1): 113-116.

    王若, 底青云, 王妙月等. 2009. 用積分方程法研究源與勘探區(qū)之間的三維體對CSAMT觀測曲線的影響. 地球物理學(xué)報, 52(6): 1573-1582, doi: 10.3969/j.issn.0001-5733.2009.06.019.

    吳桂桔, 胡祥云, 劉慧. 2010. CSAMT三維正演數(shù)值模擬研究進(jìn)展. 地球物理進(jìn)展, 25(5): 1795-1801, doi: 10.3969/j.issn.1004-2903.2010.05.037.

    徐志鋒, 吳小平. 2010. 可控源電磁三維頻率域有限元模擬. 地球物理學(xué)報, 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019.

    楊波, 徐義賢, 何展翔等. 2012. 考慮海底地形的三維頻率域可控源電磁響應(yīng)有限體積法模擬. 地球物理學(xué)報, 55(4): 1390-1399, doi: 10.6038/j.issn.0001-5733.2012.04.035.

    (本文編輯 胡素芳)

    3-D frequency-domain CSEM modeling using a parallel direct solver

    HAN Bo, HU Xiang-Yun*, HUANG Yi-Fan, PENG Rong-Hua, LI Jian-Hui, CAI Jian-Chao

    HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China

    Three-dimensional modeling of electromagnetic data is a computationally demanding problem. For frequently-used numerical techniques such as finite-element and finite-difference methods, solving the large linear systems arising from the discretization of Maxwell's equations is a key step which has a major impact on the applicability of the solution, and it has always been a research topic to solve the linear equations efficiently, robustly and accurately. A 3D modeling scheme based on direct solutions of the linear system is presented for frequency domain controlled-source electromagnetic (CSEM) surveys. The Helmholtz equation in terms of secondary electric fields is discretized using a finite-volume (FV) method over a staggered grid. Taking advantage of recent developments in numerical algorithms and the availability of computational resources, the resulting linear system of FV equations is solved directly using the massively parallel solver, namely MUMPS, instead of the most commonly used linear solvers, i.e. Krylov subspace iterative techniques. The direct solver carries out an LU (and possibly LDLT) decomposition of the system matrix and then computes solutions efficiently by applying forward and backward substitutions.To evaluate the computational performance of the direct solver, a series of numerical tests based on synthetic 1D models were conducted, and the results indicate that (1) Normalized residuals of solutions are almost independent of the conductivity value assigned to air layers but increase rapidly as the frequency value decreases. Nevertheless, the order of magnitude of the largest normalized residual is as small as 10-11. At the same time, although the matrix factorization time varies as either the air conductivity or the frequency changes, the variation is only a fraction of the total run time. (2) Both the execution time and required memory increase rapidly (more than linearly) with increasing grid sizes. (3) By executing MUMPS in parallel over multiple processors, not only the total run time but also the average memory used per processor can be reduced a lot. However, the total memory requirement increases with the number of processes. The scalability of MUMPS is limited. Additional numerical experiments considering specific survey settings were done to demonstrate the reliability and effectiveness of the code, and the results are as follows: (1) The FV numerical solutions show excellent agreement with semi-analytic solutions for the 1D models. (2) The computation time of a multitransmitter problem is comparable to that of a single-transmitter problem. (3) Reasonable modeling results of 3D models can be obtained for both typical land and marine survey scenarios.In summary, compared with iterative linear solvers, the direct solvers generally benefit CSEM modeling in three aspects. The first is they often provide more accurate solutions. The second is that direct solvers are much more stable for ill-conditioned linear systems, which are almost inevitable because of large electrical conductivity contrasts and/or non-uniform grid. The last is in multitransmitter problems, only a single matrix factorization is necessary, and multiple solutions can be achieved very easily by reusing the factors. The presented 3D CSEM modeling scheme, which employs the MUMPS direct solver, possesses all these advantages. In addition, solving linear systems can be executed in parallel to speed up the computation and to reduce the average memory used per node, although the parallel scalability of MUMPS is limited. In spite of the fact that matrix factorizations for large models can entail tremendous computational cost, it can be anticipated that direct solvers will be used more and more widely as the development of both numerical algorithms and computers.

    Controlled-source electromagnetics; 3-D modeling; LU factorization; Direct solver

    國家自然科學(xué)基金項目(41274077、41474055)和中國地質(zhì)調(diào)查局項目(12120113101800)聯(lián)合資助.

    韓波,男,1987年生,博士研究生,研究方向為電磁法的正演與反演模擬.E-mail:hanbo1735@163.com

    *通訊作者胡祥云,男,1966年生,教授,主要從事電磁法的理論與應(yīng)用研究.E-mail:xyhu@cug.edu.cn

    10.6038/cjg20150816.

    10.6038/cjg20150816

    P631

    2015-02-12,2015-06-16收修定稿

    韓波,胡祥云,黃一凡等. 2015. 基于并行化直接解法的頻率域可控源電磁三維正演.地球物理學(xué)報,58(8):2812-2826,

    Han B, Hu X Y, Huang Y F, et al. 2015. 3-D frequency-domain CSEM modeling using a parallel direct solver.ChineseJ.Geophys. (in Chinese),58(8):2812-2826,doi:10.6038/cjg20150816.

    猜你喜歡
    場源內(nèi)存電磁
    例談求解疊加電場的電場強(qiáng)度的策略
    基于深度展開ISTA網(wǎng)絡(luò)的混合源定位方法
    信號處理(2022年10期)2022-11-16 00:50:56
    基于矩陣差分的遠(yuǎn)場和近場混合源定位方法
    “春夏秋冬”的內(nèi)存
    三維多孔電磁復(fù)合支架構(gòu)建與理化表征
    掌握基礎(chǔ)知識 不懼電磁偏轉(zhuǎn)
    一種識別位場場源的混合小波方法
    電磁換向閥應(yīng)用探討
    河南科技(2014年16期)2014-02-27 14:13:21
    瞬變電磁法在煤礦防治水中的應(yīng)用
    河南科技(2014年8期)2014-02-27 14:07:41
    基于內(nèi)存的地理信息訪問技術(shù)
    深夜精品福利| 亚洲精品乱久久久久久| 久久99热这里只频精品6学生| 一级a爱视频在线免费观看| 女人被躁到高潮嗷嗷叫费观| 日韩欧美一区二区三区在线观看 | 老司机靠b影院| 考比视频在线观看| 国产有黄有色有爽视频| 日韩制服骚丝袜av| 少妇猛男粗大的猛烈进出视频| 午夜两性在线视频| 精品乱码久久久久久99久播| 免费看十八禁软件| 一本色道久久久久久精品综合| 国产成人免费无遮挡视频| 丰满饥渴人妻一区二区三| 国产亚洲av片在线观看秒播厂| 女人久久www免费人成看片| 亚洲精品自拍成人| 国产日韩一区二区三区精品不卡| 韩国高清视频一区二区三区| 高清视频免费观看一区二区| 一区二区三区精品91| 国产精品 欧美亚洲| 免费日韩欧美在线观看| 天天影视国产精品| 纵有疾风起免费观看全集完整版| 午夜成年电影在线免费观看| 精品少妇一区二区三区视频日本电影| av视频免费观看在线观看| 狠狠狠狠99中文字幕| 中文字幕人妻熟女乱码| 亚洲中文日韩欧美视频| 丝袜喷水一区| 99国产精品一区二区蜜桃av | 波多野结衣一区麻豆| 国产高清videossex| 久久久久网色| 日韩精品免费视频一区二区三区| 美女大奶头黄色视频| 十分钟在线观看高清视频www| 中文字幕高清在线视频| 精品亚洲成国产av| 国产免费福利视频在线观看| 成人三级做爰电影| 欧美国产精品va在线观看不卡| www.自偷自拍.com| 99精国产麻豆久久婷婷| 2018国产大陆天天弄谢| 国产精品免费视频内射| 欧美午夜高清在线| tocl精华| 精品国产国语对白av| 嫩草影视91久久| 19禁男女啪啪无遮挡网站| 欧美精品一区二区免费开放| 一级毛片电影观看| a在线观看视频网站| 久久精品久久久久久噜噜老黄| 国产男人的电影天堂91| 99香蕉大伊视频| 亚洲av电影在线进入| 欧美精品一区二区大全| 最黄视频免费看| 自拍欧美九色日韩亚洲蝌蚪91| 久久ye,这里只有精品| 丰满人妻熟妇乱又伦精品不卡| 精品欧美一区二区三区在线| 国产成人啪精品午夜网站| 人人妻人人添人人爽欧美一区卜| 99久久综合免费| 老熟妇乱子伦视频在线观看 | 国产一卡二卡三卡精品| 成人免费观看视频高清| 亚洲av日韩在线播放| tocl精华| 免费久久久久久久精品成人欧美视频| 久久青草综合色| 欧美人与性动交α欧美软件| 建设人人有责人人尽责人人享有的| 国产深夜福利视频在线观看| kizo精华| 男女免费视频国产| 亚洲人成电影免费在线| 成人影院久久| 国产精品久久久久久精品古装| www.999成人在线观看| 欧美日韩视频精品一区| 一区二区三区乱码不卡18| videos熟女内射| 国内毛片毛片毛片毛片毛片| 黑人猛操日本美女一级片| 九色亚洲精品在线播放| 美女中出高潮动态图| 性高湖久久久久久久久免费观看| 亚洲激情五月婷婷啪啪| 亚洲成国产人片在线观看| 丁香六月欧美| 一本久久精品| 动漫黄色视频在线观看| 午夜影院在线不卡| 精品人妻熟女毛片av久久网站| 性色av乱码一区二区三区2| 亚洲一卡2卡3卡4卡5卡精品中文| 亚洲国产欧美在线一区| 久久99热这里只频精品6学生| 一边摸一边抽搐一进一出视频| 亚洲全国av大片| 久热这里只有精品99| 精品国产乱码久久久久久男人| 美女主播在线视频| 亚洲欧洲精品一区二区精品久久久| 亚洲专区字幕在线| 亚洲欧美色中文字幕在线| 欧美在线一区亚洲| 亚洲欧美日韩高清在线视频 | 亚洲精品美女久久av网站| 精品熟女少妇八av免费久了| 不卡av一区二区三区| 日韩电影二区| 如日韩欧美国产精品一区二区三区| 国产av国产精品国产| 欧美+亚洲+日韩+国产| 99精国产麻豆久久婷婷| 成人18禁高潮啪啪吃奶动态图| 国产无遮挡羞羞视频在线观看| 在线天堂中文资源库| 日本wwww免费看| 夜夜骑夜夜射夜夜干| 男女无遮挡免费网站观看| 亚洲一区二区三区欧美精品| 丁香六月欧美| 日韩视频在线欧美| 在线观看免费午夜福利视频| 超碰成人久久| 日本精品一区二区三区蜜桃| 在线av久久热| 亚洲中文字幕日韩| 国产av精品麻豆| 91麻豆av在线| 亚洲国产欧美日韩在线播放| 亚洲精品中文字幕一二三四区 | 性色av一级| 在线天堂中文资源库| 久久久久久亚洲精品国产蜜桃av| 久久国产精品影院| 亚洲av日韩精品久久久久久密| 女人精品久久久久毛片| 妹子高潮喷水视频| 亚洲精品中文字幕一二三四区 | 在线观看www视频免费| 国产三级黄色录像| 亚洲avbb在线观看| 亚洲欧美激情在线| 黄色视频,在线免费观看| 亚洲国产欧美在线一区| 飞空精品影院首页| 高清黄色对白视频在线免费看| 国产亚洲午夜精品一区二区久久| 亚洲专区中文字幕在线| 黄频高清免费视频| 中国国产av一级| 欧美+亚洲+日韩+国产| 亚洲国产精品一区二区三区在线| 高潮久久久久久久久久久不卡| 老熟女久久久| 亚洲av国产av综合av卡| 亚洲精品国产av成人精品| 热99久久久久精品小说推荐| 日韩制服骚丝袜av| 丰满饥渴人妻一区二区三| 国产精品.久久久| 成人国产av品久久久| 亚洲专区中文字幕在线| 美女大奶头黄色视频| 黑人欧美特级aaaaaa片| 嫩草影视91久久| 99国产精品一区二区蜜桃av | 夫妻午夜视频| 亚洲伊人色综图| 狠狠婷婷综合久久久久久88av| 欧美国产精品一级二级三级| 国产野战对白在线观看| 香蕉国产在线看| 亚洲av成人一区二区三| 他把我摸到了高潮在线观看 | 91精品伊人久久大香线蕉| 欧美精品高潮呻吟av久久| 这个男人来自地球电影免费观看| 欧美乱码精品一区二区三区| 久久中文看片网| 三级毛片av免费| 色老头精品视频在线观看| 午夜福利,免费看| 亚洲国产av影院在线观看| 久久久久网色| 90打野战视频偷拍视频| 国产成人啪精品午夜网站| 国产欧美日韩一区二区三区在线| 最新的欧美精品一区二区| 国产成人a∨麻豆精品| 国产精品香港三级国产av潘金莲| 黄色片一级片一级黄色片| 国产老妇伦熟女老妇高清| 人妻一区二区av| 性少妇av在线| 亚洲精品国产区一区二| 最近最新中文字幕大全免费视频| 黄片播放在线免费| 国产精品香港三级国产av潘金莲| av片东京热男人的天堂| 欧美亚洲日本最大视频资源| 伊人亚洲综合成人网| 亚洲色图综合在线观看| 午夜福利一区二区在线看| 亚洲欧美成人综合另类久久久| 欧美日韩精品网址| 黑丝袜美女国产一区| 欧美另类一区| 国产亚洲精品第一综合不卡| 后天国语完整版免费观看| 久久精品亚洲熟妇少妇任你| 国产欧美亚洲国产| 亚洲avbb在线观看| 精品亚洲乱码少妇综合久久| 两性夫妻黄色片| 亚洲欧美清纯卡通| 成年人免费黄色播放视频| 国产精品免费大片| 91字幕亚洲| 久久精品成人免费网站| 搡老熟女国产l中国老女人| 国产精品国产av在线观看| 国产麻豆69| av免费在线观看网站| 黑丝袜美女国产一区| 久久热在线av| 在线观看一区二区三区激情| 三级毛片av免费| 99热全是精品| 丝袜美腿诱惑在线| 91老司机精品| 考比视频在线观看| 18禁裸乳无遮挡动漫免费视频| 自拍欧美九色日韩亚洲蝌蚪91| 国产精品亚洲av一区麻豆| 超碰97精品在线观看| 人成视频在线观看免费观看| 脱女人内裤的视频| 午夜福利视频在线观看免费| 国产成人精品久久二区二区91| 中文欧美无线码| 一本—道久久a久久精品蜜桃钙片| 国产成人精品久久二区二区免费| 免费在线观看影片大全网站| 久久中文看片网| 亚洲五月婷婷丁香| 成人影院久久| 国产精品久久久久久精品古装| 成年人免费黄色播放视频| 国产国语露脸激情在线看| 在线观看免费高清a一片| 成年av动漫网址| 丝袜美腿诱惑在线| 人妻久久中文字幕网| 国产亚洲av高清不卡| 亚洲久久久国产精品| 香蕉国产在线看| 久久影院123| 国产日韩欧美视频二区| 欧美黑人欧美精品刺激| 亚洲美女黄色视频免费看| 窝窝影院91人妻| 丝袜在线中文字幕| av超薄肉色丝袜交足视频| 日韩免费高清中文字幕av| 91精品三级在线观看| 777米奇影视久久| 18禁观看日本| tocl精华| 欧美日韩黄片免| 国产精品 国内视频| 国产成+人综合+亚洲专区| 老汉色av国产亚洲站长工具| 制服人妻中文乱码| 成年动漫av网址| 首页视频小说图片口味搜索| 久久精品国产a三级三级三级| tocl精华| 亚洲精品中文字幕在线视频| 久久九九热精品免费| 亚洲国产精品成人久久小说| 精品免费久久久久久久清纯 | 国产精品免费大片| tube8黄色片| 久久这里只有精品19| 啦啦啦啦在线视频资源| 国产成人影院久久av| 久久久久网色| 亚洲av日韩在线播放| 精品久久久精品久久久| 俄罗斯特黄特色一大片| 亚洲中文日韩欧美视频| 91成年电影在线观看| 黄片播放在线免费| 99精品欧美一区二区三区四区| 2018国产大陆天天弄谢| 成年动漫av网址| 成年女人毛片免费观看观看9 | 丝袜美足系列| 亚洲中文av在线| 色视频在线一区二区三区| 午夜免费观看性视频| 亚洲成国产人片在线观看| 精品亚洲乱码少妇综合久久| 性少妇av在线| 国产野战对白在线观看| 肉色欧美久久久久久久蜜桃| 搡老熟女国产l中国老女人| 国产一区二区激情短视频 | 国产在线免费精品| 啦啦啦在线免费观看视频4| 国产不卡av网站在线观看| 人妻久久中文字幕网| 欧美av亚洲av综合av国产av| 精品久久蜜臀av无| 国产真人三级小视频在线观看| 国产成人免费观看mmmm| 午夜福利视频在线观看免费| 老司机午夜十八禁免费视频| 亚洲伊人色综图| 一进一出抽搐动态| 日韩制服丝袜自拍偷拍| 日韩制服骚丝袜av| av在线app专区| 欧美日韩一级在线毛片| 亚洲全国av大片| av福利片在线| 韩国精品一区二区三区| 国产熟女午夜一区二区三区| 免费在线观看影片大全网站| 中文字幕人妻熟女乱码| 自线自在国产av| 久9热在线精品视频| 热re99久久精品国产66热6| 亚洲激情五月婷婷啪啪| 午夜两性在线视频| 亚洲精华国产精华精| 久久久久网色| 亚洲精华国产精华精| 水蜜桃什么品种好| 大型av网站在线播放| 一区二区三区乱码不卡18| 久久国产精品大桥未久av| 精品卡一卡二卡四卡免费| 岛国毛片在线播放| 免费一级毛片在线播放高清视频 | 亚洲欧美成人综合另类久久久| 老熟妇仑乱视频hdxx| 国产一区二区三区av在线| 久久人人爽人人片av| 亚洲精品自拍成人| 天堂俺去俺来也www色官网| 亚洲激情五月婷婷啪啪| 一区二区三区激情视频| a在线观看视频网站| 欧美成人午夜精品| 日本av免费视频播放| 婷婷色av中文字幕| 99热国产这里只有精品6| 飞空精品影院首页| 一区二区av电影网| 人人妻,人人澡人人爽秒播| av福利片在线| 999久久久精品免费观看国产| 中国美女看黄片| 热99国产精品久久久久久7| 人人妻人人爽人人添夜夜欢视频| 久久中文字幕一级| 夜夜骑夜夜射夜夜干| 黄频高清免费视频| 成人影院久久| 91麻豆精品激情在线观看国产 | 80岁老熟妇乱子伦牲交| 人妻久久中文字幕网| 大片电影免费在线观看免费| 国产日韩欧美亚洲二区| 99精品久久久久人妻精品| 午夜激情av网站| 黄频高清免费视频| 亚洲国产欧美网| 国产在线视频一区二区| 成人18禁高潮啪啪吃奶动态图| 日本欧美视频一区| 最黄视频免费看| 亚洲中文字幕日韩| 久久天躁狠狠躁夜夜2o2o| 超碰成人久久| 丝袜美腿诱惑在线| 真人做人爱边吃奶动态| 91麻豆精品激情在线观看国产 | 日韩一区二区三区影片| 人人妻人人澡人人爽人人夜夜| 国产精品一区二区在线不卡| av天堂在线播放| 最近最新免费中文字幕在线| 老熟妇乱子伦视频在线观看 | 老司机福利观看| 成人av一区二区三区在线看 | 亚洲三区欧美一区| 精品免费久久久久久久清纯 | 国内毛片毛片毛片毛片毛片| 国产精品一区二区免费欧美 | 爱豆传媒免费全集在线观看| 国产无遮挡羞羞视频在线观看| 亚洲欧美一区二区三区黑人| 午夜福利视频精品| 精品乱码久久久久久99久播| 2018国产大陆天天弄谢| 十八禁高潮呻吟视频| 国产一区二区三区av在线| 欧美精品亚洲一区二区| 亚洲精品中文字幕一二三四区 | 曰老女人黄片| 成在线人永久免费视频| 91精品国产国语对白视频| 国产精品.久久久| 搡老熟女国产l中国老女人| 亚洲精品中文字幕一二三四区 | 丰满饥渴人妻一区二区三| 日本av免费视频播放| 亚洲精品国产av蜜桃| 国产精品秋霞免费鲁丝片| 午夜福利视频在线观看免费| 免费av中文字幕在线| 1024香蕉在线观看| 青青草视频在线视频观看| 久久国产精品男人的天堂亚洲| 国产av又大| www日本在线高清视频| 国内毛片毛片毛片毛片毛片| 啦啦啦在线免费观看视频4| 欧美97在线视频| 国产极品粉嫩免费观看在线| 欧美性长视频在线观看| 久久久久精品人妻al黑| www.熟女人妻精品国产| av在线老鸭窝| 亚洲一卡2卡3卡4卡5卡精品中文| 丝袜脚勾引网站| 国产伦人伦偷精品视频| 美女福利国产在线| 男女床上黄色一级片免费看| 伦理电影免费视频| 在线观看免费视频网站a站| 国产成人免费观看mmmm| www.自偷自拍.com| 最黄视频免费看| 婷婷色av中文字幕| 欧美亚洲 丝袜 人妻 在线| 久9热在线精品视频| 欧美亚洲日本最大视频资源| e午夜精品久久久久久久| 日韩精品免费视频一区二区三区| 天堂8中文在线网| 亚洲av电影在线进入| 成年人午夜在线观看视频| 桃红色精品国产亚洲av| av线在线观看网站| 成人三级做爰电影| 亚洲久久久国产精品| 久久精品国产综合久久久| 欧美在线一区亚洲| 美女视频免费永久观看网站| 首页视频小说图片口味搜索| 国产不卡av网站在线观看| 啦啦啦视频在线资源免费观看| 最近中文字幕2019免费版| 欧美激情久久久久久爽电影 | 正在播放国产对白刺激| 欧美日韩中文字幕国产精品一区二区三区 | 电影成人av| 国产一级毛片在线| 少妇的丰满在线观看| 人妻 亚洲 视频| av在线播放精品| 性色av乱码一区二区三区2| 国产视频一区二区在线看| 一区在线观看完整版| 午夜福利视频在线观看免费| 成人18禁高潮啪啪吃奶动态图| 久久久国产精品麻豆| 成年人午夜在线观看视频| 亚洲美女黄色视频免费看| 新久久久久国产一级毛片| 久久久精品94久久精品| www.av在线官网国产| 夜夜骑夜夜射夜夜干| 亚洲精品久久久久久婷婷小说| 成年av动漫网址| 久久久国产一区二区| 在线精品无人区一区二区三| 免费在线观看视频国产中文字幕亚洲 | 一区二区三区乱码不卡18| 麻豆av在线久日| 91字幕亚洲| 亚洲avbb在线观看| 欧美激情极品国产一区二区三区| av不卡在线播放| 两性夫妻黄色片| 在线观看人妻少妇| 搡老乐熟女国产| 成人免费观看视频高清| 国产精品一区二区精品视频观看| 一进一出抽搐动态| 日韩欧美一区视频在线观看| 在线观看免费高清a一片| 亚洲黑人精品在线| 动漫黄色视频在线观看| 各种免费的搞黄视频| 国产又色又爽无遮挡免| 国产成人a∨麻豆精品| av网站免费在线观看视频| 纯流量卡能插随身wifi吗| 高清av免费在线| 91精品三级在线观看| 超碰成人久久| 国产欧美日韩一区二区精品| 精品国产一区二区三区四区第35| 亚洲色图综合在线观看| 亚洲av欧美aⅴ国产| 天天躁夜夜躁狠狠躁躁| 美女扒开内裤让男人捅视频| 在线天堂中文资源库| 亚洲一区中文字幕在线| 欧美日韩福利视频一区二区| 天天影视国产精品| 啦啦啦中文免费视频观看日本| 久久99一区二区三区| 国产高清视频在线播放一区 | 男女床上黄色一级片免费看| av不卡在线播放| 麻豆国产av国片精品| 国产一区有黄有色的免费视频| 国产成人av激情在线播放| 亚洲七黄色美女视频| 婷婷色av中文字幕| 老司机福利观看| 国产在线免费精品| 亚洲欧美精品综合一区二区三区| 亚洲国产av影院在线观看| 国产又色又爽无遮挡免| 一本久久精品| 欧美日韩视频精品一区| 国产淫语在线视频| 亚洲中文日韩欧美视频| 汤姆久久久久久久影院中文字幕| av在线app专区| 国产精品一区二区免费欧美 | 国产无遮挡羞羞视频在线观看| 在线av久久热| 国产在线视频一区二区| 欧美日韩精品网址| 极品人妻少妇av视频| 精品乱码久久久久久99久播| 亚洲va日本ⅴa欧美va伊人久久 | 日韩一区二区三区影片| 各种免费的搞黄视频| 成年人免费黄色播放视频| 777久久人妻少妇嫩草av网站| 一进一出抽搐动态| 国产不卡av网站在线观看| 亚洲欧美清纯卡通| 69精品国产乱码久久久| 久久久精品免费免费高清| 咕卡用的链子| 欧美日韩视频精品一区| 在线观看免费视频网站a站| 王馨瑶露胸无遮挡在线观看| 亚洲人成电影观看| 免费黄频网站在线观看国产| 亚洲色图综合在线观看| 亚洲一区二区三区欧美精品| 一个人免费看片子| 国精品久久久久久国模美| 男女边摸边吃奶| 午夜福利视频在线观看免费| av在线播放精品| 视频在线观看一区二区三区| tube8黄色片| 啦啦啦视频在线资源免费观看| 亚洲av国产av综合av卡| 亚洲欧美色中文字幕在线| 久久 成人 亚洲| 精品福利永久在线观看| 黑人巨大精品欧美一区二区蜜桃| 久久热在线av| 人成视频在线观看免费观看| av天堂在线播放| 黄色怎么调成土黄色| 69av精品久久久久久 | 久久人妻福利社区极品人妻图片| 热re99久久精品国产66热6| 免费久久久久久久精品成人欧美视频| 久久国产精品大桥未久av| 婷婷成人精品国产| 丰满人妻熟妇乱又伦精品不卡| 少妇裸体淫交视频免费看高清 | 男人操女人黄网站| av欧美777| 一本综合久久免费| 一区在线观看完整版| 另类亚洲欧美激情| 又黄又粗又硬又大视频|