• <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ù)
    一级毛片女人18水好多| 国产精品日韩av在线免费观看| 精品久久久久久,| 动漫黄色视频在线观看| 亚洲精品一卡2卡三卡4卡5卡| 又黄又粗又硬又大视频| 一区二区三区国产精品乱码| 国产片内射在线| tocl精华| 国产精品久久视频播放| 黄色视频,在线免费观看| 久久精品aⅴ一区二区三区四区| 女人高潮潮喷娇喘18禁视频| 亚洲人成网站在线播放欧美日韩| 又黄又粗又硬又大视频| 久久久久国产精品人妻aⅴ院| 男男h啪啪无遮挡| 99精品在免费线老司机午夜| 淫妇啪啪啪对白视频| 精品国产乱子伦一区二区三区| 日韩欧美国产在线观看| 正在播放国产对白刺激| 亚洲avbb在线观看| 国产v大片淫在线免费观看| 欧美成人免费av一区二区三区| 99热这里只有精品一区 | 一进一出抽搐gif免费好疼| 宅男免费午夜| 国产av又大| 午夜久久久久精精品| 777久久人妻少妇嫩草av网站| 婷婷精品国产亚洲av在线| 国产一区二区激情短视频| 给我免费播放毛片高清在线观看| 国产成人av激情在线播放| 欧美成人午夜精品| 天堂动漫精品| 可以在线观看的亚洲视频| 欧美一级a爱片免费观看看 | 亚洲av美国av| 午夜福利视频1000在线观看| 精品国产乱子伦一区二区三区| 妹子高潮喷水视频| 亚洲男人天堂网一区| 视频在线观看一区二区三区| 欧美久久黑人一区二区| 国产午夜福利久久久久久| 久久天堂一区二区三区四区| 欧美成狂野欧美在线观看| 国产高清videossex| 777久久人妻少妇嫩草av网站| 久久亚洲精品不卡| 一个人观看的视频www高清免费观看 | 亚洲午夜精品一区,二区,三区| 波多野结衣av一区二区av| 90打野战视频偷拍视频| 国产成人一区二区三区免费视频网站| 一本综合久久免费| 国产三级在线视频| 日韩av在线大香蕉| svipshipincom国产片| 成年版毛片免费区| 久久中文看片网| 久久久久国内视频| 中文字幕久久专区| 免费在线观看影片大全网站| 国产精品一区二区三区四区久久 | 久久中文字幕人妻熟女| 日本熟妇午夜| 午夜激情av网站| 国产精品国产高清国产av| 亚洲va日本ⅴa欧美va伊人久久| 亚洲成人免费电影在线观看| 日韩三级视频一区二区三区| 99精品久久久久人妻精品| 高潮久久久久久久久久久不卡| 国产成人欧美| 国产精品自产拍在线观看55亚洲| 首页视频小说图片口味搜索| 亚洲va日本ⅴa欧美va伊人久久| 欧美在线黄色| 亚洲国产精品成人综合色| 美女扒开内裤让男人捅视频| 日本成人三级电影网站| 亚洲精品美女久久av网站| 亚洲一卡2卡3卡4卡5卡精品中文| 中文资源天堂在线| 日韩 欧美 亚洲 中文字幕| 女性生殖器流出的白浆| 亚洲激情在线av| 丝袜在线中文字幕| 色综合婷婷激情| 国产一级毛片七仙女欲春2 | 久久热在线av| 午夜福利在线在线| 国产精品爽爽va在线观看网站 | 国产片内射在线| 在线观看66精品国产| 国产免费av片在线观看野外av| 欧美zozozo另类| 熟女少妇亚洲综合色aaa.| 超碰成人久久| 一区二区日韩欧美中文字幕| 国产成人影院久久av| 欧美国产日韩亚洲一区| 国产一区二区三区视频了| 精品乱码久久久久久99久播| 亚洲一卡2卡3卡4卡5卡精品中文| 国产精品野战在线观看| 美女免费视频网站| 亚洲五月婷婷丁香| 成人一区二区视频在线观看| 日本五十路高清| 亚洲狠狠婷婷综合久久图片| 精品久久蜜臀av无| 激情在线观看视频在线高清| 91麻豆精品激情在线观看国产| 一个人观看的视频www高清免费观看 | 亚洲av成人一区二区三| 神马国产精品三级电影在线观看 | 男女视频在线观看网站免费 | 亚洲精品一卡2卡三卡4卡5卡| 免费在线观看完整版高清| 精品熟女少妇八av免费久了| 在线国产一区二区在线| 久久99热这里只有精品18| 男人的好看免费观看在线视频 | 午夜视频精品福利| 露出奶头的视频| 可以免费在线观看a视频的电影网站| 老司机在亚洲福利影院| 国产一卡二卡三卡精品| 好男人在线观看高清免费视频 | av视频在线观看入口| 级片在线观看| 伊人久久大香线蕉亚洲五| 丝袜美腿诱惑在线| 国产av一区在线观看免费| 亚洲精品在线美女| 国内毛片毛片毛片毛片毛片| 日韩欧美在线二视频| 国产在线精品亚洲第一网站| 欧美+亚洲+日韩+国产| 国产一区二区三区在线臀色熟女| 久久久久久亚洲精品国产蜜桃av| 国产亚洲精品久久久久久毛片| 亚洲av成人av| 欧美一级毛片孕妇| 色哟哟哟哟哟哟| 久久久久久久午夜电影| 在线播放国产精品三级| 999久久久国产精品视频| 一夜夜www| 日本在线视频免费播放| 亚洲精品在线观看二区| 99国产精品一区二区蜜桃av| 黄色毛片三级朝国网站| 国产主播在线观看一区二区| 午夜激情福利司机影院| av免费在线观看网站| 久久精品影院6| 亚洲免费av在线视频| 国产精品乱码一区二三区的特点| 十八禁网站免费在线| 很黄的视频免费| 色婷婷久久久亚洲欧美| 国产av不卡久久| 国产区一区二久久| 日本黄色视频三级网站网址| 制服丝袜大香蕉在线| 国产97色在线日韩免费| 黄色成人免费大全| 老汉色av国产亚洲站长工具| 久久久久国产一级毛片高清牌| 国产精品一区二区精品视频观看| 成人亚洲精品一区在线观看| 中文字幕人妻丝袜一区二区| 国产成人影院久久av| 久久 成人 亚洲| 亚洲国产毛片av蜜桃av| tocl精华| 国产精品99久久99久久久不卡| 99在线人妻在线中文字幕| 97超级碰碰碰精品色视频在线观看| 美女国产高潮福利片在线看| 美女高潮到喷水免费观看| 老司机福利观看| 色哟哟哟哟哟哟| 黄频高清免费视频| 可以在线观看毛片的网站| 国产精品久久久av美女十八| 午夜免费激情av| 在线国产一区二区在线| 一进一出好大好爽视频| 亚洲一区高清亚洲精品| 不卡一级毛片| 性色av乱码一区二区三区2| 两个人看的免费小视频| 啦啦啦免费观看视频1| 十八禁人妻一区二区| 亚洲国产日韩欧美精品在线观看 | 人妻丰满熟妇av一区二区三区| 日韩精品青青久久久久久| 亚洲av熟女| 啦啦啦韩国在线观看视频| 亚洲欧美日韩无卡精品| 在线观看日韩欧美| 国产精品乱码一区二三区的特点| 国产午夜精品久久久久久| 免费看日本二区| 欧美三级亚洲精品| 欧美zozozo另类| 国产精品电影一区二区三区| 国产精品久久久人人做人人爽| 亚洲成人久久爱视频| 欧美国产精品va在线观看不卡| 老汉色av国产亚洲站长工具| 亚洲天堂国产精品一区在线| 国产又色又爽无遮挡免费看| 国产熟女午夜一区二区三区| 亚洲专区国产一区二区| 日本a在线网址| 黄色a级毛片大全视频| 国产又爽黄色视频| 热99re8久久精品国产| 少妇 在线观看| 午夜福利免费观看在线| 免费在线观看日本一区| av欧美777| 久久久精品国产亚洲av高清涩受| 亚洲精品在线美女| 色精品久久人妻99蜜桃| 又黄又粗又硬又大视频| 淫秽高清视频在线观看| 我的亚洲天堂| 国产蜜桃级精品一区二区三区| 99久久国产精品久久久| 妹子高潮喷水视频| 日韩大码丰满熟妇| 欧美日韩精品网址| 欧美激情高清一区二区三区| 日韩大尺度精品在线看网址| 一进一出好大好爽视频| 色综合欧美亚洲国产小说| 不卡一级毛片| 午夜免费激情av| 天天一区二区日本电影三级| 欧美日韩中文字幕国产精品一区二区三区| 国内精品久久久久久久电影| 国产在线精品亚洲第一网站| 超碰成人久久| 国产成人系列免费观看| 69av精品久久久久久| 精品久久久久久久人妻蜜臀av| 久久婷婷人人爽人人干人人爱| 精品无人区乱码1区二区| 国产精品99久久99久久久不卡| av中文乱码字幕在线| 啦啦啦观看免费观看视频高清| 中文字幕高清在线视频| 国产一级毛片七仙女欲春2 | 国产成人av教育| 美国免费a级毛片| 中文字幕人成人乱码亚洲影| 中亚洲国语对白在线视频| 91大片在线观看| 夜夜看夜夜爽夜夜摸| 国产av不卡久久| 1024视频免费在线观看| 啪啪无遮挡十八禁网站| 久久伊人香网站| 日韩欧美国产一区二区入口| 在线观看午夜福利视频| 免费观看精品视频网站| 丁香欧美五月| 午夜福利视频1000在线观看| 精品卡一卡二卡四卡免费| 一区二区日韩欧美中文字幕| 久久香蕉激情| 搞女人的毛片| 日韩欧美三级三区| 久久精品91无色码中文字幕| 亚洲精品av麻豆狂野| 波多野结衣av一区二区av| 国产精品 国内视频| 成人精品一区二区免费| 天天躁狠狠躁夜夜躁狠狠躁| 免费人成视频x8x8入口观看| 欧美日韩亚洲国产一区二区在线观看| 麻豆成人av在线观看| 午夜免费观看网址| 午夜免费激情av| 欧美最黄视频在线播放免费| ponron亚洲| 欧美日韩亚洲综合一区二区三区_| 校园春色视频在线观看| 久久香蕉精品热| 国产亚洲欧美在线一区二区| e午夜精品久久久久久久| 少妇熟女aⅴ在线视频| 在线永久观看黄色视频| 国产色视频综合| 香蕉久久夜色| 村上凉子中文字幕在线| 国产成人av教育| 身体一侧抽搐| 麻豆国产av国片精品| 亚洲人成网站在线播放欧美日韩| 这个男人来自地球电影免费观看| 亚洲五月色婷婷综合| 国产高清有码在线观看视频 | 久久国产亚洲av麻豆专区| 亚洲国产欧洲综合997久久, | 美女高潮到喷水免费观看| xxx96com| 亚洲激情在线av| 欧洲精品卡2卡3卡4卡5卡区| av在线播放免费不卡| 国产成人欧美在线观看| 美女大奶头视频| 国产97色在线日韩免费| 精品久久久久久,| 国产日本99.免费观看| 91字幕亚洲| 亚洲一卡2卡3卡4卡5卡精品中文| 国产99白浆流出| 俺也久久电影网| 又紧又爽又黄一区二区| a在线观看视频网站| 亚洲av五月六月丁香网| 亚洲精华国产精华精| 嫩草影院精品99| 国内久久婷婷六月综合欲色啪| 国产精品久久久av美女十八| 中出人妻视频一区二区| 日韩 欧美 亚洲 中文字幕| 国产精品亚洲av一区麻豆| www.自偷自拍.com| aaaaa片日本免费| 国内少妇人妻偷人精品xxx网站 | 日韩高清综合在线| 又紧又爽又黄一区二区| 国产精品 国内视频| 欧美不卡视频在线免费观看 | 色综合欧美亚洲国产小说| 亚洲成人精品中文字幕电影| 精品一区二区三区四区五区乱码| 亚洲av成人不卡在线观看播放网| 午夜激情福利司机影院| 2021天堂中文幕一二区在线观 | 久久精品国产亚洲av香蕉五月| 欧美一级毛片孕妇| 久久精品夜夜夜夜夜久久蜜豆 | 国产成人av激情在线播放| 国产精品,欧美在线| 亚洲国产精品久久男人天堂| 久久天堂一区二区三区四区| 久久久久久久久免费视频了| 日本精品一区二区三区蜜桃| 成人欧美大片| 亚洲狠狠婷婷综合久久图片| 又黄又粗又硬又大视频| 欧美zozozo另类| 激情在线观看视频在线高清| 欧美成人一区二区免费高清观看 | 超碰成人久久| 日韩av在线大香蕉| www日本在线高清视频| 免费观看人在逋| 女人爽到高潮嗷嗷叫在线视频| 欧美性猛交╳xxx乱大交人| 韩国av一区二区三区四区| 90打野战视频偷拍视频| 国产精品亚洲美女久久久| 日韩欧美免费精品| 两性午夜刺激爽爽歪歪视频在线观看 | 久久中文字幕一级| 色婷婷久久久亚洲欧美| 色精品久久人妻99蜜桃| 精品国产亚洲在线| 亚洲人成电影免费在线| 国产av在哪里看| 在线看三级毛片| 国产三级在线视频| 国产亚洲欧美精品永久| 草草在线视频免费看| 高清在线国产一区| 亚洲成人国产一区在线观看| 久久青草综合色| 久久久水蜜桃国产精品网| 日日摸夜夜添夜夜添小说| 19禁男女啪啪无遮挡网站| 啦啦啦韩国在线观看视频| 美女高潮到喷水免费观看| 在线国产一区二区在线| 搡老妇女老女人老熟妇| 国产精品1区2区在线观看.| 久久久久国内视频| 老司机靠b影院| 18禁黄网站禁片午夜丰满| 国产精品影院久久| 久久久国产成人精品二区| 日韩成人在线观看一区二区三区| 国产精品久久电影中文字幕| 亚洲欧美日韩无卡精品| 国产精品美女特级片免费视频播放器 | 久久亚洲精品不卡| 最新在线观看一区二区三区| 日韩免费av在线播放| 夜夜爽天天搞| 欧美国产精品va在线观看不卡| 日本精品一区二区三区蜜桃| 成年免费大片在线观看| 日韩欧美 国产精品| 亚洲五月婷婷丁香| 曰老女人黄片| 少妇裸体淫交视频免费看高清 | 国产久久久一区二区三区| 神马国产精品三级电影在线观看 | 丝袜人妻中文字幕| 成人欧美大片| 亚洲国产看品久久| 欧美日韩福利视频一区二区| 欧美日韩精品网址| 国产一卡二卡三卡精品| 三级毛片av免费| 国产精品1区2区在线观看.| 国内毛片毛片毛片毛片毛片| 欧美最黄视频在线播放免费| 国产免费av片在线观看野外av| 97超级碰碰碰精品色视频在线观看| 妹子高潮喷水视频| 日日干狠狠操夜夜爽| 久久久久久久久中文| 天堂动漫精品| 午夜久久久在线观看| 精品久久久久久,| 日韩有码中文字幕| 亚洲精品在线美女| 国产精品久久视频播放| 日韩国内少妇激情av| 国产一区二区激情短视频| 人人澡人人妻人| 亚洲片人在线观看| 法律面前人人平等表现在哪些方面| 在线永久观看黄色视频| 啦啦啦观看免费观看视频高清| 午夜视频精品福利| 亚洲一卡2卡3卡4卡5卡精品中文| 午夜福利在线在线| 亚洲av电影不卡..在线观看| 国产私拍福利视频在线观看| 成人一区二区视频在线观看| 制服丝袜大香蕉在线| 欧美不卡视频在线免费观看 | 中文在线观看免费www的网站 | 别揉我奶头~嗯~啊~动态视频| 午夜福利在线在线| 法律面前人人平等表现在哪些方面| 国产私拍福利视频在线观看| 精品高清国产在线一区| 亚洲av第一区精品v没综合| 欧美精品亚洲一区二区| 给我免费播放毛片高清在线观看| 伊人久久大香线蕉亚洲五| 狂野欧美激情性xxxx| av在线播放免费不卡| 两个人看的免费小视频| 成人亚洲精品一区在线观看| 美女国产高潮福利片在线看| 午夜精品在线福利| av天堂在线播放| 亚洲国产中文字幕在线视频| 日本精品一区二区三区蜜桃| 国产亚洲精品一区二区www| 久久久久久久精品吃奶| 丝袜在线中文字幕| 国产97色在线日韩免费| 亚洲激情在线av| 真人做人爱边吃奶动态| 久久精品国产亚洲av高清一级| 亚洲三区欧美一区| 日韩有码中文字幕| av在线播放免费不卡| 级片在线观看| 午夜福利免费观看在线| 精品福利观看| 国产欧美日韩一区二区精品| 亚洲国产精品999在线| 中文在线观看免费www的网站 | 国产日本99.免费观看| 亚洲第一av免费看| 夜夜躁狠狠躁天天躁| 亚洲九九香蕉| 亚洲成人免费电影在线观看| aaaaa片日本免费| 午夜免费鲁丝| 精品国产乱码久久久久久男人| 亚洲av五月六月丁香网| 在线观看日韩欧美| 无遮挡黄片免费观看| 国产精品久久电影中文字幕| 国产精品二区激情视频| 熟妇人妻久久中文字幕3abv| 欧美成人一区二区免费高清观看 | 亚洲一区二区三区色噜噜| 亚洲欧美激情综合另类| 精品国产美女av久久久久小说| 我的亚洲天堂| 叶爱在线成人免费视频播放| 国产精品98久久久久久宅男小说| 可以在线观看的亚洲视频| 美女午夜性视频免费| 精品国产乱子伦一区二区三区| 少妇裸体淫交视频免费看高清 | 亚洲成人久久性| 首页视频小说图片口味搜索| 18禁国产床啪视频网站| 亚洲av片天天在线观看| 丝袜人妻中文字幕| 免费看a级黄色片| 午夜久久久久精精品| 国产一区二区三区在线臀色熟女| 婷婷精品国产亚洲av| 一级作爱视频免费观看| 国产在线观看jvid| 日本一本二区三区精品| 欧美日韩亚洲综合一区二区三区_| 黑人巨大精品欧美一区二区mp4| 国产国语露脸激情在线看| 国产亚洲精品第一综合不卡| xxx96com| 欧美日韩亚洲综合一区二区三区_| 熟女电影av网| 在线av久久热| 在线免费观看的www视频| 精品高清国产在线一区| 欧美成人一区二区免费高清观看 | 999精品在线视频| aaaaa片日本免费| 国产91精品成人一区二区三区| 中文资源天堂在线| 嫁个100分男人电影在线观看| ponron亚洲| 午夜成年电影在线免费观看| 天堂动漫精品| 亚洲午夜精品一区,二区,三区| 制服诱惑二区| 美女免费视频网站| 亚洲在线自拍视频| 黄频高清免费视频| 欧美三级亚洲精品| 欧美一级毛片孕妇| 亚洲久久久国产精品| 欧美日韩瑟瑟在线播放| 国产亚洲av嫩草精品影院| 波多野结衣av一区二区av| 一本精品99久久精品77| 听说在线观看完整版免费高清| 激情在线观看视频在线高清| 国产成人av教育| 久久亚洲精品不卡| 精品国产国语对白av| 午夜福利在线观看吧| 欧美又色又爽又黄视频| 欧美另类亚洲清纯唯美| 国产精品九九99| 婷婷亚洲欧美| 亚洲人成伊人成综合网2020| 日韩欧美一区二区三区在线观看| 好看av亚洲va欧美ⅴa在| 久久久久国内视频| 亚洲欧美日韩无卡精品| 成人一区二区视频在线观看| 久久亚洲真实| 国产男靠女视频免费网站| 天堂影院成人在线观看| 精品久久久久久久毛片微露脸| 日韩中文字幕欧美一区二区| 很黄的视频免费| 久久久久久久久久黄片| 又黄又粗又硬又大视频| 两性夫妻黄色片| 国产一区在线观看成人免费| 啦啦啦韩国在线观看视频| 精品久久久久久久末码| 色在线成人网| 国产又黄又爽又无遮挡在线| 久久精品国产亚洲av高清一级| 午夜视频精品福利| 国产高清videossex| 黑丝袜美女国产一区| 一区二区三区激情视频| 十分钟在线观看高清视频www| 久久精品影院6| 国产精品国产高清国产av| x7x7x7水蜜桃| 午夜福利一区二区在线看| 又黄又粗又硬又大视频| 午夜精品在线福利| 一区二区日韩欧美中文字幕| av电影中文网址| 精品熟女少妇八av免费久了| 91麻豆精品激情在线观看国产| 无限看片的www在线观看| 可以在线观看的亚洲视频| 宅男免费午夜| 美女国产高潮福利片在线看| 日本撒尿小便嘘嘘汇集6| 精品国产美女av久久久久小说| 日韩国内少妇激情av| 可以在线观看的亚洲视频| 欧美午夜高清在线| 欧美绝顶高潮抽搐喷水| 国产成人av教育|