MT有限元模擬中截?cái)噙吔绲挠绊?/h1>
2013-09-22 13:50:30湯井田
湯井田，薛 帥
中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083
0 引言
國內(nèi)外學(xué)者很早就進(jìn)行了大地電磁場(chǎng)的數(shù)值模擬研究。1971年Coggon［1］首先提出將有限單元法應(yīng)用在電磁法數(shù)值模擬中，并詳細(xì)闡述了模擬大地電磁場(chǎng)的有限元格式。有限單元法作為一種便捷通用的數(shù)值分析方法，相較于有限差分法和積分方程法，具有易于擬合復(fù)雜結(jié)構(gòu)、易于對(duì)復(fù)雜介質(zhì)目標(biāo)建模、程序通用性強(qiáng)等優(yōu)勢(shì)，從而被迅速廣泛地應(yīng)用于電磁場(chǎng)的數(shù)值模擬中。國外學(xué)者Rijo［2］、Wannamaker等［3－4］分別對(duì)電磁場(chǎng)的有限元模擬進(jìn)行了發(fā)展和完善，Wannamaker等［3－4］將有限元中的矩形網(wǎng)格繼續(xù)剖分為4個(gè)三角形網(wǎng)格，可以更好地模擬地形和地下構(gòu)造，提高了數(shù)值計(jì)算精度。在國內(nèi)，1981年，陳樂壽［4］首先詳細(xì)介紹了有限元在大地電磁場(chǎng)正演計(jì)算中的應(yīng)用，并利用矩形－三角形網(wǎng)格剖分對(duì)有限元模擬大地電磁進(jìn)行了改進(jìn)。之后胡建德等［6］、徐世浙等［7］、楊長福［8］、陳小斌［9］分別對(duì)地球物理場(chǎng)中的有限單元法進(jìn)行了研究和發(fā)展，至今，有限單元法已經(jīng)成為二維和三維大地電磁數(shù)值模擬的主要手段。
在利用有限元方法模擬二維和三維大地電磁場(chǎng)時(shí)，電磁場(chǎng)正演計(jì)算時(shí)的一些邊界條件需要在“無窮遠(yuǎn)處”才能夠滿足。而有限元的模擬區(qū)域是有限的，所以網(wǎng)格邊界處存在截?cái)噙吔缬绊憽Ｔ谟邢拊M大地電磁場(chǎng)時(shí)，為了避免截?cái)噙吔绲挠绊懀话銓⒕W(wǎng)格邊界放在距離研究區(qū)域外很遠(yuǎn)處，如 Rijo［2］、Wannamaker［3－4］、陳樂壽等［10］、徐世浙［11］、趙廣茂等［12］、劉長生等［13］在進(jìn)行大地電磁有限元數(shù)值模擬時(shí)網(wǎng)格邊界都放在足夠遠(yuǎn)處來滿足邊界條件。但是“足夠遠(yuǎn)距離”一直沒有適當(dāng)?shù)膮⒖贾担W(wǎng)格邊界選擇小了會(huì)影響計(jì)算精度，而選擇過大則需要較高的內(nèi)存和付出較長的計(jì)算時(shí)間。Sharma等［14］曾利用趨膚深度δ來進(jìn)行網(wǎng)格剖分和劃分網(wǎng)格邊界，認(rèn)為將網(wǎng)格邊界設(shè)置在5δ處時(shí)即可忽略截?cái)噙吔绲挠绊懀菦]有進(jìn)行相應(yīng)的研究說明。作者利用有限元二維大地電磁正演程序，計(jì)算總結(jié)出適合有限元大地電磁模擬的參考網(wǎng)格邊界。
1 邊界條件和截?cái)噙吔?/h2>
假定地下電性結(jié)構(gòu)為二維結(jié)構(gòu)，選取右旋直角坐標(biāo)系的原點(diǎn)在地面上，z軸正方向垂直向下，x軸的方向平行于走向方向，即介質(zhì)的電性參數(shù)僅與y、z方向有關(guān)。根據(jù)麥克斯韋方程，二維地電結(jié)構(gòu)下的大地電磁滿足以下二式，稱之為亥娒霍茲方程式［10］：
其中：Hx和Ex是x 方向磁場(chǎng)和電場(chǎng)；k2＝iωμσ；ω是角頻率；μ是導(dǎo)磁率；σ是導(dǎo)電率。
圖1 E型波（a）和 H型波（b）的研究區(qū)域［11］Fig.1 Research area of E type wave（a）and H type wave（b）［11］
邊界如圖1所示，設(shè)u為需要求解的電磁場(chǎng)，在上邊界AB處，電磁場(chǎng)滿足第一類邊界條件；值得注意的是，TE（transverse electric）模式的上邊界要離開地面一段距離。在左右邊界AD和BC處，要求由橫向不均勻體引起的二次場(chǎng)為0。為此，有2種方式的邊界條件加載［10］：一種是強(qiáng)加邊界條件，即計(jì)算出水平層狀介質(zhì)中的電磁場(chǎng)分布來給出邊界上的值；另一種是采用齊次的第二類邊界條件，即在側(cè)邊界上的場(chǎng)分量u的法向?qū)?shù)為零。一般采用后一種方式作為側(cè)邊界的邊界條件。而下邊界CD處以下被認(rèn)為是均質(zhì)巖石，要求局部異常體在CD上的異常場(chǎng)為0。所以u(píng)滿足的邊界條件為
當(dāng)?shù)叵麓嬖诓痪鶆蝮w時(shí)，網(wǎng)格邊界需要在足夠遠(yuǎn)處才能夠滿足電磁場(chǎng)的邊界條件。
2 截?cái)噙吔缪芯?/h2>
由于在利用大地電磁（magnetotellurics，MT）有限元正演程序進(jìn)行大地電磁數(shù)值模擬時(shí)，網(wǎng)格剖分會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生很大的影響，所以在進(jìn)行截?cái)噙吔缪芯繒r(shí)，要先避免因網(wǎng)格剖分而造成的計(jì)算結(jié)果差異。筆者為確保計(jì)算結(jié)果的正確與可靠，在橫向上（y方向）－1 000～1 000m、縱向上（z方向）0～2 000m，按照10m的網(wǎng)格間距進(jìn)行剖分；而在該研究區(qū)域以外，無論橫向還是縱向的網(wǎng)格間距，皆按1.1倍率進(jìn)行增長剖分。對(duì)于TE模式參照Wannamaker等［3－4］在PW2D中的作法，以總網(wǎng)格橫向距離的一半作為上網(wǎng)格邊界的距離。
數(shù)值模擬時(shí)，在不改變網(wǎng)格剖分的情況下，只對(duì)外網(wǎng)格邊界進(jìn)行擴(kuò)大和減小，從而避免網(wǎng)格剖分的影響。根據(jù)電磁場(chǎng)理論，某一頻率的電磁場(chǎng)在該頻點(diǎn)的一個(gè)趨膚深度內(nèi)迅速衰減至1／e，并在傳播中繼續(xù)衰減至0；因此，筆者以某一頻點(diǎn)的趨膚深度作為該頻點(diǎn)網(wǎng)格邊界選擇的量度：
式中：δ為趨膚深度；ρ為電阻率；f為頻率?？梢?，δ與f相關(guān)。為了體現(xiàn)普遍性，在MT儀器常用的4個(gè)頻段中分別選擇一個(gè)頻率作為研究頻率，即0.25、2.50、25.00、250.00Hz。
在模型選擇上，筆者分別對(duì)一維模型（均勻半空間模型和層狀介質(zhì)模型）和二維地電模型進(jìn)行模擬計(jì)算。在進(jìn)行數(shù)值模擬前，將2個(gè)方向的網(wǎng)格邊界分別放在20倍的趨膚深度距離處，計(jì)算出一組數(shù)據(jù)作為參考解，認(rèn)為該參考解是在“足夠遠(yuǎn)”的網(wǎng)格邊界條件下計(jì)算的較精確的數(shù)值解。
2.1 一維模型模擬
首先對(duì)電阻率為100Ω·m的均勻半空間模型進(jìn)行模擬計(jì)算，分別計(jì)算左右網(wǎng)格邊界和下網(wǎng)格邊界同時(shí)為100m和同時(shí)為500m的情況，計(jì)算結(jié)果如表1。
再建立一個(gè)層狀介質(zhì)模型（圖2），分別計(jì)算左右網(wǎng)格邊界和下網(wǎng)格邊界同時(shí)為500m和同時(shí)為1 000m情況下的視電阻率，結(jié)果如表2所示。
表1 均勻半空間視電阻率計(jì)算結(jié)果Table1 Apparent resistivity results of homogenous half space
表2 H型層狀介質(zhì)視電阻率計(jì)算結(jié)果Table2 Apparent resistivity results of layered medium model
在一維模型中不存在局部異常體，自然滿足有限元模擬的外邊界條件。通過對(duì)以上2個(gè)模型的計(jì)算，進(jìn)一步說明了在一維有限元數(shù)值模擬時(shí)邊界條件自然滿足，不存在截?cái)噙吔缬绊?，從而無需進(jìn)行網(wǎng)格邊界的擴(kuò)展。
圖2 層狀介質(zhì)模型Fig.2 Layered medium model
2.2 二維地電模型
有限元模擬地下存在二度體的二維地電模型時(shí)，由于局部不均勻體的存在，不合適的網(wǎng)格邊界會(huì)產(chǎn)生較大的誤差。對(duì)模型2D－1（圖3）進(jìn)行模擬計(jì)算，設(shè)置左右網(wǎng)格邊界和下網(wǎng)格邊界均為500m，分別在O點(diǎn)和B點(diǎn)進(jìn)行觀測(cè)，表3是O點(diǎn)觀測(cè)結(jié)果。
從表3中可以看出，不合適的網(wǎng)格邊界對(duì)計(jì)算結(jié)果影響較大。因此，對(duì)二維地電模型進(jìn)行有限元模擬時(shí)，需要選擇合適的網(wǎng)格邊界。
圖3 二維地電模型Fig.3 Two－dimensional geoeletric model
根據(jù)所選頻率和模型電阻率，將下網(wǎng)格邊界設(shè)置在“足夠遠(yuǎn)”的距離，放在200 000m的地方，計(jì)算趨膚深度作為變化量度來改變左右網(wǎng)格邊界的邊界距離，將計(jì)算結(jié)果與參考解進(jìn)行比較。同時(shí)約定下網(wǎng)格邊界以地下異常體的下邊界為基準(zhǔn)進(jìn)行擴(kuò)展，左右網(wǎng)格邊界以異常體的左右邊界為基準(zhǔn)進(jìn)行擴(kuò)展。約定計(jì)算趨膚深度時(shí)的電阻率為地下電性結(jié)構(gòu)中電阻率最高者。
對(duì)背景介質(zhì)中存在低阻異常體和高阻異常體2種模型進(jìn)行模擬計(jì)算，并在二度體中心上方O點(diǎn)和邊界上方B點(diǎn)同時(shí)進(jìn)行TE模式和TM（transverse magnetic）模式觀測(cè)。
2.2.1 左右網(wǎng)格邊界變化
表3 模型2D－1在O點(diǎn)觀測(cè)結(jié)果Table3 Observed results of model 2D－1on point O
表4 左右網(wǎng)格邊界變化下視電阻率觀測(cè)結(jié)果Table4 Apparent resistivity results when changing left and right mesh boundary
從表4可看出：對(duì)于模型2D－1，TE模式在1δ網(wǎng)格邊界時(shí)，O點(diǎn)和B點(diǎn)觀測(cè)結(jié)果的相對(duì)誤差除250Hz外大部降至1%以下；在2δ的網(wǎng)格邊界時(shí)，O點(diǎn)觀測(cè)結(jié)果的相對(duì)誤差已全降至0.1%以下，B點(diǎn)雖相對(duì)于1δ網(wǎng)格邊界條件時(shí)的相對(duì)誤差有所減小，但部分頻率的誤差仍在0.1%數(shù)量級(jí)上；繼續(xù)增大網(wǎng)格邊界至3δ時(shí)，O點(diǎn)和B點(diǎn)觀測(cè)的相對(duì)誤差已全部降至0.1%以下；繼續(xù)加大網(wǎng)格邊界至5δ和10δ，觀測(cè)結(jié)果與參考結(jié)果相等，誤差為0。
對(duì)于TM模式：在1δ網(wǎng)格邊界時(shí)，計(jì)算的結(jié)果相對(duì)誤差全部在1%以下；增大到2δ網(wǎng)格邊界時(shí)，相對(duì)誤差大部為0；繼續(xù)加大網(wǎng)格邊界，相對(duì)誤差保持為0。對(duì)于模型2D－2，在1δ網(wǎng)格邊界時(shí)，大部分計(jì)算結(jié)果與參考解相等，相對(duì)誤差為0。
值得注意的是，對(duì)模型2D－2進(jìn)行模擬時(shí)，網(wǎng)格邊界變化的量度采用的趨膚深度已相當(dāng)于模型2D－1的3倍。
2.2.2 下網(wǎng)格邊界變化
在對(duì)下邊界研究時(shí)，則把左右網(wǎng)格邊界設(shè)置在“足夠遠(yuǎn)”的距離，以趨膚深度作為量度來改變下網(wǎng)格邊界的邊界距離，計(jì)算結(jié)果并與參考解進(jìn)行比較。仍然以模型2D－1和模型2D－2為計(jì)算模型。模型2D－1的部分計(jì)算結(jié)果如表5。
從表5中看出：根據(jù)TE模式計(jì)算結(jié)果，在1δ網(wǎng)格邊界時(shí)，O點(diǎn)和B點(diǎn)觀測(cè)結(jié)果的相對(duì)誤差已全部降至1%以下；在2δ網(wǎng)格邊界時(shí)，相對(duì)誤差繼續(xù)減小，但大部仍在0.1%的數(shù)量級(jí)上；繼續(xù)增大網(wǎng)格邊界至3δ時(shí)，O點(diǎn)和B點(diǎn)觀測(cè)的相對(duì)誤差已全部降至0.1%以下；加大邊界至5δ和10δ，觀測(cè)結(jié)果與參考結(jié)果相等，誤差為0。
對(duì)于TM模式：在1δ網(wǎng)格邊界時(shí)，計(jì)算的結(jié)果相對(duì)誤差大部在1%以下；增大到2δ網(wǎng)格邊界時(shí)，除B點(diǎn)觀測(cè)的25Hz和250Hz外，大部分結(jié)果的相對(duì)誤差降至0.1%以下；繼續(xù)增大網(wǎng)格邊界，誤差繼續(xù)減少直至為0。
對(duì)于下網(wǎng)格邊界變化下的高阻模型2D－2，計(jì)算的結(jié)果與對(duì)應(yīng)的左右網(wǎng)格邊界的變化結(jié)果相似。即在1δ時(shí)，大部分計(jì)算結(jié)果已與參考解相等，相對(duì)誤差為0。
3 截?cái)噙吔邕x擇
3.1 綜合網(wǎng)格邊界計(jì)算
以上的計(jì)算結(jié)論都是在某一方向的網(wǎng)格邊界位于“足夠遠(yuǎn)處”得到的，現(xiàn)在把以上得出的最佳網(wǎng)格邊界距離綜合在一起進(jìn)行計(jì)算比較。
對(duì)于模型2D－1：計(jì)算TE模式時(shí)，左右網(wǎng)格邊界和下網(wǎng)格邊界都放在3δ處；計(jì)算TM模式時(shí)，則將左右網(wǎng)格邊界和下網(wǎng)格邊界放在2δ處。對(duì)于模型2D－2，把網(wǎng)格邊界全部設(shè)置在1δ處即可。綜合網(wǎng)格邊界條件下部分觀測(cè)結(jié)果如表6所示。
由表6可見，以上利用綜合網(wǎng)格邊界計(jì)算的結(jié)果相對(duì)誤差基本在0.1%以下，說明該網(wǎng)格邊界條件下計(jì)算結(jié)果的可靠性。
3.2 綜合模型計(jì)算
為了進(jìn)一步驗(yàn)證綜合的最佳網(wǎng)格邊界計(jì)算結(jié)果的可靠性，建立1個(gè)較復(fù)雜的二維模型，并采用新的頻率進(jìn)行計(jì)算。模型2D－3（圖4）是在背景介質(zhì)中賦存3種異常體的較復(fù)雜模型，分別在O、O1、O2和B點(diǎn)觀測(cè)，計(jì)算的頻率為0.01、0.10、1.00、10.00、100.00、1 000.00Hz，剖分情況與前面相同。
在模型2D－3中存在多個(gè)異常體，在確定網(wǎng)格邊界變化量度時(shí)，需要幾個(gè)邊界變化量度進(jìn)行比較才能確定。在TE模式計(jì)算時(shí)，對(duì)于低阻體ρ2和ρ4，當(dāng)它們賦存在較高電阻率的背景介質(zhì)中時(shí)，趨膚深度計(jì)算采用背景電阻率，邊界變化量度為3δ＝503；而當(dāng)背景介質(zhì)中存在高阻體ρ3時(shí)，趨膚深度計(jì)算采用高阻體電阻率ρ3，其邊界變化量度為δ ＝503-。為了確保計(jì)算結(jié)果正確和保證精度，網(wǎng)格邊界變化量度應(yīng)為TM模式計(jì)算時(shí)網(wǎng)格邊界變化量度的確定與之類似。表7為模型2D－3的部分計(jì)算結(jié)果。
表5 模型2D－1下網(wǎng)格邊界變化下O點(diǎn)視電阻率觀測(cè)結(jié)果Table5 Apparent resistivity results of model 2D－1on point Owhen changing bottom boundary
表6 綜合網(wǎng)格邊界O點(diǎn)觀測(cè)結(jié)果Table6 Observed results on point Owith compositive boundaries
表7 模型2D－3在O點(diǎn)的觀測(cè)結(jié)果Table7 Observed results of model 2D－3on point O
圖4 模型2D－3Fig.4 Model 2D－3
由表7可見，利用綜合的網(wǎng)格邊界計(jì)算的結(jié)果與參考解的相對(duì)誤差全部在0.1%以下，進(jìn)一步驗(yàn)證了綜合邊界的正確性。
值得注意的是，以上的結(jié)果是在背景電阻率和異常體電阻率最大相差一個(gè)數(shù)量級(jí)的情況下得出的。作者又利用綜合邊界計(jì)算了背景電阻率為100 Ω·m、異常體電阻率分別為1、0.1、0.01、0.001 Ω·m的情況以及背景電阻率為100Ω·m、異常體電阻率為200Ω·m的情況，計(jì)算結(jié)果表明TE模式中誤差最大為0.40%，TM模式誤差基本為0。另外作者還計(jì)算了一些極端的情況，如背景電阻率為1Ω·m、異常體電阻率為2Ω·m的模型和背景電阻率為2Ω·m、異常體電阻率為1Ω·m的模型，發(fā)現(xiàn)參考邊界在此也適用。
4 結(jié)論與建議
利用有限元方法對(duì)大地電磁進(jìn)行模擬計(jì)算時(shí)，由于存在截?cái)噙吔缬绊懀瑥亩绊懹?jì)算結(jié)果和精度，或占用大量內(nèi)存花費(fèi)較長的時(shí)間。筆者通過大量的數(shù)值模擬計(jì)算與對(duì)比，以趨膚深度作為網(wǎng)格邊界變化量度，考慮到大地電磁的精度要求，總結(jié)得出適合二維大地電磁有限元正演模擬計(jì)算的參考網(wǎng)格邊界：
1）利用有限元對(duì)一維地電模型進(jìn)行模擬計(jì)算時(shí)，不存在截?cái)噙吔绲挠绊?，邊界條件自然滿足。
2）對(duì)二維低阻模型進(jìn)行模擬計(jì)算時(shí)，以背景電阻率計(jì)算趨膚深度作為網(wǎng)格邊界變化量度。對(duì)于TE模式，當(dāng)左右網(wǎng)格邊界和下網(wǎng)格邊界都放在3δ處時(shí)，截?cái)噙吔绲挠绊懸芽珊雎浴?duì)于TM模式，左右網(wǎng)格邊界和下網(wǎng)格邊界在2δ處，截?cái)噙吔缬绊懣梢院雎浴?/p>
3）對(duì)二維高阻模型進(jìn)行模擬計(jì)算時(shí)，以高阻體電阻率計(jì)算趨膚深度作為網(wǎng)格邊界變化量度。不管TE模式還是TM模式，當(dāng)左右網(wǎng)格邊界和下網(wǎng)格邊界位于1δ時(shí)，截?cái)噙吔缬绊懸芽珊雎浴?/p>
4）當(dāng)二維地電模型中，同時(shí)存在低阻異常體和高阻異常體時(shí)，分別以背景介質(zhì)電阻率和高阻體中電阻率最高者計(jì)算趨膚深度，然后根據(jù)參考網(wǎng)格邊界選擇最佳的網(wǎng)格邊界。
（References）：
［1］Coggon J H.Electromagnetic and Electrical Modeling by Finite Element Method［J］.Geophysics，1971，36（1）：132－155.
［2］Rodi W L.A Technique for Improving the Accuracy of Finite Element Solution for Magnetotelluric Data［J］.Geophysics，1976，44（2）：483－506.
［3］Wannamaker P E，Stodt J A，Rijo W L.Two－Dimensional Topographic Responses in Magnetotelluric Modeling Using Finite Element［J］.Geophysics，1986，51（11）：2131－2144.
［4］Wannamaker P E，Stodt J A，Rijo W L.A Stable Finite Element Solution for Two－Dimensional Magnetotelluric Modelling［J］.Geophysics，1987，88（1）：277－296.
［5］陳樂壽.有限元法在大地電磁場(chǎng)正演計(jì)算中的應(yīng)用與改進(jìn)［J］.石油物探，1981，20（3）：84－103.Chen Leshou.The Improvement and Application of Finite Element Method to 2DForward Solution of Magnetotelluric Sounding［J］.Geophysical Prospecting for Petroleum，1981，20（3）：84－103.
［6］胡建德，蔡剛.用三角形二次插值有限元法計(jì)算二維大地電磁測(cè)深曲線［J］.石油地球物理勘探，1984，19（4）：358－367.Hu Jiande，Cai Gang.Calculating 2DMagnetotelluric Sounding Curve with the Twice Interpolation Triangle Finite Element Method ［J］. Oil Geophysical Prospecting，1984，19（4）：358－367.
［7］史明娟，徐世浙，劉斌.大地電磁二次函數(shù)插值的有限元法正演模擬［J］.地球物理學(xué)報(bào)，1997，40（3）：421－430.Shi Mingjuan，Xu Shizhe，Liu Bin.Finite Element Method Using Quadratic Element in MT Forward Modeling［J］.Acta Geophysica Sinica，1997，40（3）：421－430.
［8］楊長福.大地電磁二維對(duì)稱各向異性介質(zhì)的有限元數(shù)值模擬［J］.西北地震學(xué)報(bào)，1997，19（2）：27－33.Yang Changfu.MT Numerical Simulation of Symmetrically 2DAnisotropic Media Based on the Finite Element Method ［J］. Northwestern Seismological Journal，1997，19（2）：27－33.
［9］陳小斌.有限元直接迭代算法［J］.物探化探計(jì)算技術(shù)，1999，21（2）：165－171.Chen Xiaobin.Direct Iterative Algorithm of Finite Element［J］.Computing Techniques for Geophysical and Geochemical Exploration，1999，21（2）：165－171.
［10］陳樂壽，劉任，王天生.大地電磁測(cè)深資料處理與解釋［M］.北京：石油工業(yè)出版社，1989：1－160.Chen Leshou， Liu Ren， Wang Tiansheng.Magnetotelluric Sounding Data Processing and Interpretation［M］.Beijing：Petroleum Industry Press，1989：1－160.
［11］徐世浙.地球物理中的有限單元法［M］.北京：科學(xué)出版社，1994：220－241.Xu Shizhe.The Finite Element in Geophysics［M］.Beijing：Science Press，1994：220－241.
［12］趙廣茂，李桐林，王大勇，等.基于二次場(chǎng)二維起伏地形MT有限元數(shù)值模擬［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2008，38（6）：1055－1059.Zhao Guangmao，Li Tonglin，Wang Dayong，et al.Secondary Field－Based Two－Dimensional Topographic Numerical Simulation in Magnetoteilurics by Finite Element Method［J］.Journal of Jilin University：Earth Science Edition，2008，38（6）：1055－1059.
［13］劉長生，湯井田，任政勇，等.基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬［J］.中南大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，41（5）：1855－1859.Liu Changsheng，Tang Jingtian，Ren Zhengyong，et al.Three－Dimension Magnetotellurics Modeling by Adaptive Edge Finite－Element Using Unstructured Meshes［J］.Journal of Central South University：Science and Technology，2010，41（5）：1855－1859.
［14］Sharma P S，Kaikkonen P.An Automated Finite Element Mesh Generation and Element Coding in 2－D Electromagnetic Inversion［J］.Geophysics，1998，34（3）：93－114.
猜你喜歡
一類帶有Stieltjes積分邊界條件的分?jǐn)?shù)階微分方程邊值問題正解應(yīng)用數(shù)學(xué)(2020年4期)2020-12-28 00:36:38 帶有積分邊界條件的奇異攝動(dòng)邊值問題的漸近解數(shù)學(xué)物理學(xué)報(bào)(2020年5期)2020-11-26 06:06:30 不等高軟橫跨橫向承力索計(jì)算及計(jì)算結(jié)果判斷研究甘肅科技(2020年20期)2020-04-13 00:30:40 三維電阻率成像與高聚物注漿在水閘加固中的應(yīng)用水利科技與經(jīng)濟(jì)(2016年3期)2016-04-22 01:05:00 隨鉆電阻率測(cè)井的固定探測(cè)深度合成方法西南石油大學(xué)學(xué)報(bào)(自然科學(xué)版)(2015年4期)2015-08-20 09:05:16 海洋可控源電磁場(chǎng)視電阻率計(jì)算方法物探化探計(jì)算技術(shù)(2015年2期)2015-02-28 17:42:42 帶Robin邊界條件的2維隨機(jī)Ginzburg-Landau方程的吸引子四川師范大學(xué)學(xué)報(bào)（自然科學(xué)版）(2015年1期)2015-02-28 14:07:21 超壓測(cè)試方法對(duì)炸藥TNT當(dāng)量計(jì)算結(jié)果的影響火炸藥學(xué)報(bào)(2014年3期)2014-03-20 13:17:39 粉煤灰摻量對(duì)水泥漿體電阻率與自收縮的影響建筑材料學(xué)報(bào)(2014年3期)2014-03-11 17:08:01 帶非齊次邊界條件的p—Laplacian方程正解的存在唯一性東北師大學(xué)報(bào)（自然科學(xué)版）(2014年1期)2014-02-27 08:02:09
湯井田，薛 帥
中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083
0 引言
國內(nèi)外學(xué)者很早就進(jìn)行了大地電磁場(chǎng)的數(shù)值模擬研究。1971年Coggon［1］首先提出將有限單元法應(yīng)用在電磁法數(shù)值模擬中，并詳細(xì)闡述了模擬大地電磁場(chǎng)的有限元格式。有限單元法作為一種便捷通用的數(shù)值分析方法，相較于有限差分法和積分方程法，具有易于擬合復(fù)雜結(jié)構(gòu)、易于對(duì)復(fù)雜介質(zhì)目標(biāo)建模、程序通用性強(qiáng)等優(yōu)勢(shì)，從而被迅速廣泛地應(yīng)用于電磁場(chǎng)的數(shù)值模擬中。國外學(xué)者Rijo［2］、Wannamaker等［3－4］分別對(duì)電磁場(chǎng)的有限元模擬進(jìn)行了發(fā)展和完善，Wannamaker等［3－4］將有限元中的矩形網(wǎng)格繼續(xù)剖分為4個(gè)三角形網(wǎng)格，可以更好地模擬地形和地下構(gòu)造，提高了數(shù)值計(jì)算精度。在國內(nèi)，1981年，陳樂壽［4］首先詳細(xì)介紹了有限元在大地電磁場(chǎng)正演計(jì)算中的應(yīng)用，并利用矩形－三角形網(wǎng)格剖分對(duì)有限元模擬大地電磁進(jìn)行了改進(jìn)。之后胡建德等［6］、徐世浙等［7］、楊長福［8］、陳小斌［9］分別對(duì)地球物理場(chǎng)中的有限單元法進(jìn)行了研究和發(fā)展，至今，有限單元法已經(jīng)成為二維和三維大地電磁數(shù)值模擬的主要手段。
在利用有限元方法模擬二維和三維大地電磁場(chǎng)時(shí)，電磁場(chǎng)正演計(jì)算時(shí)的一些邊界條件需要在“無窮遠(yuǎn)處”才能夠滿足。而有限元的模擬區(qū)域是有限的，所以網(wǎng)格邊界處存在截?cái)噙吔缬绊憽Ｔ谟邢拊M大地電磁場(chǎng)時(shí)，為了避免截?cái)噙吔绲挠绊懀话銓⒕W(wǎng)格邊界放在距離研究區(qū)域外很遠(yuǎn)處，如 Rijo［2］、Wannamaker［3－4］、陳樂壽等［10］、徐世浙［11］、趙廣茂等［12］、劉長生等［13］在進(jìn)行大地電磁有限元數(shù)值模擬時(shí)網(wǎng)格邊界都放在足夠遠(yuǎn)處來滿足邊界條件。但是“足夠遠(yuǎn)距離”一直沒有適當(dāng)?shù)膮⒖贾担W(wǎng)格邊界選擇小了會(huì)影響計(jì)算精度，而選擇過大則需要較高的內(nèi)存和付出較長的計(jì)算時(shí)間。Sharma等［14］曾利用趨膚深度δ來進(jìn)行網(wǎng)格剖分和劃分網(wǎng)格邊界，認(rèn)為將網(wǎng)格邊界設(shè)置在5δ處時(shí)即可忽略截?cái)噙吔绲挠绊懀菦]有進(jìn)行相應(yīng)的研究說明。作者利用有限元二維大地電磁正演程序，計(jì)算總結(jié)出適合有限元大地電磁模擬的參考網(wǎng)格邊界。
1 邊界條件和截?cái)噙吔?/h2>
假定地下電性結(jié)構(gòu)為二維結(jié)構(gòu)，選取右旋直角坐標(biāo)系的原點(diǎn)在地面上，z軸正方向垂直向下，x軸的方向平行于走向方向，即介質(zhì)的電性參數(shù)僅與y、z方向有關(guān)。根據(jù)麥克斯韋方程，二維地電結(jié)構(gòu)下的大地電磁滿足以下二式，稱之為亥娒霍茲方程式［10］：
其中：Hx和Ex是x 方向磁場(chǎng)和電場(chǎng)；k2＝iωμσ；ω是角頻率；μ是導(dǎo)磁率；σ是導(dǎo)電率。
圖1 E型波（a）和 H型波（b）的研究區(qū)域［11］Fig.1 Research area of E type wave（a）and H type wave（b）［11］
邊界如圖1所示，設(shè)u為需要求解的電磁場(chǎng)，在上邊界AB處，電磁場(chǎng)滿足第一類邊界條件；值得注意的是，TE（transverse electric）模式的上邊界要離開地面一段距離。在左右邊界AD和BC處，要求由橫向不均勻體引起的二次場(chǎng)為0。為此，有2種方式的邊界條件加載［10］：一種是強(qiáng)加邊界條件，即計(jì)算出水平層狀介質(zhì)中的電磁場(chǎng)分布來給出邊界上的值；另一種是采用齊次的第二類邊界條件，即在側(cè)邊界上的場(chǎng)分量u的法向?qū)?shù)為零。一般采用后一種方式作為側(cè)邊界的邊界條件。而下邊界CD處以下被認(rèn)為是均質(zhì)巖石，要求局部異常體在CD上的異常場(chǎng)為0。所以u(píng)滿足的邊界條件為
當(dāng)?shù)叵麓嬖诓痪鶆蝮w時(shí)，網(wǎng)格邊界需要在足夠遠(yuǎn)處才能夠滿足電磁場(chǎng)的邊界條件。
2 截?cái)噙吔缪芯?/h2>
由于在利用大地電磁（magnetotellurics，MT）有限元正演程序進(jìn)行大地電磁數(shù)值模擬時(shí)，網(wǎng)格剖分會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生很大的影響，所以在進(jìn)行截?cái)噙吔缪芯繒r(shí)，要先避免因網(wǎng)格剖分而造成的計(jì)算結(jié)果差異。筆者為確保計(jì)算結(jié)果的正確與可靠，在橫向上（y方向）－1 000～1 000m、縱向上（z方向）0～2 000m，按照10m的網(wǎng)格間距進(jìn)行剖分；而在該研究區(qū)域以外，無論橫向還是縱向的網(wǎng)格間距，皆按1.1倍率進(jìn)行增長剖分。對(duì)于TE模式參照Wannamaker等［3－4］在PW2D中的作法，以總網(wǎng)格橫向距離的一半作為上網(wǎng)格邊界的距離。
數(shù)值模擬時(shí)，在不改變網(wǎng)格剖分的情況下，只對(duì)外網(wǎng)格邊界進(jìn)行擴(kuò)大和減小，從而避免網(wǎng)格剖分的影響。根據(jù)電磁場(chǎng)理論，某一頻率的電磁場(chǎng)在該頻點(diǎn)的一個(gè)趨膚深度內(nèi)迅速衰減至1／e，并在傳播中繼續(xù)衰減至0；因此，筆者以某一頻點(diǎn)的趨膚深度作為該頻點(diǎn)網(wǎng)格邊界選擇的量度：
式中：δ為趨膚深度；ρ為電阻率；f為頻率?？梢?，δ與f相關(guān)。為了體現(xiàn)普遍性，在MT儀器常用的4個(gè)頻段中分別選擇一個(gè)頻率作為研究頻率，即0.25、2.50、25.00、250.00Hz。
在模型選擇上，筆者分別對(duì)一維模型（均勻半空間模型和層狀介質(zhì)模型）和二維地電模型進(jìn)行模擬計(jì)算。在進(jìn)行數(shù)值模擬前，將2個(gè)方向的網(wǎng)格邊界分別放在20倍的趨膚深度距離處，計(jì)算出一組數(shù)據(jù)作為參考解，認(rèn)為該參考解是在“足夠遠(yuǎn)”的網(wǎng)格邊界條件下計(jì)算的較精確的數(shù)值解。
2.1 一維模型模擬
首先對(duì)電阻率為100Ω·m的均勻半空間模型進(jìn)行模擬計(jì)算，分別計(jì)算左右網(wǎng)格邊界和下網(wǎng)格邊界同時(shí)為100m和同時(shí)為500m的情況，計(jì)算結(jié)果如表1。
再建立一個(gè)層狀介質(zhì)模型（圖2），分別計(jì)算左右網(wǎng)格邊界和下網(wǎng)格邊界同時(shí)為500m和同時(shí)為1 000m情況下的視電阻率，結(jié)果如表2所示。
表1 均勻半空間視電阻率計(jì)算結(jié)果Table1 Apparent resistivity results of homogenous half space
表2 H型層狀介質(zhì)視電阻率計(jì)算結(jié)果Table2 Apparent resistivity results of layered medium model
在一維模型中不存在局部異常體，自然滿足有限元模擬的外邊界條件。通過對(duì)以上2個(gè)模型的計(jì)算，進(jìn)一步說明了在一維有限元數(shù)值模擬時(shí)邊界條件自然滿足，不存在截?cái)噙吔缬绊?，從而無需進(jìn)行網(wǎng)格邊界的擴(kuò)展。
圖2 層狀介質(zhì)模型Fig.2 Layered medium model
2.2 二維地電模型
有限元模擬地下存在二度體的二維地電模型時(shí)，由于局部不均勻體的存在，不合適的網(wǎng)格邊界會(huì)產(chǎn)生較大的誤差。對(duì)模型2D－1（圖3）進(jìn)行模擬計(jì)算，設(shè)置左右網(wǎng)格邊界和下網(wǎng)格邊界均為500m，分別在O點(diǎn)和B點(diǎn)進(jìn)行觀測(cè)，表3是O點(diǎn)觀測(cè)結(jié)果。
從表3中可以看出，不合適的網(wǎng)格邊界對(duì)計(jì)算結(jié)果影響較大。因此，對(duì)二維地電模型進(jìn)行有限元模擬時(shí)，需要選擇合適的網(wǎng)格邊界。
圖3 二維地電模型Fig.3 Two－dimensional geoeletric model
根據(jù)所選頻率和模型電阻率，將下網(wǎng)格邊界設(shè)置在“足夠遠(yuǎn)”的距離，放在200 000m的地方，計(jì)算趨膚深度作為變化量度來改變左右網(wǎng)格邊界的邊界距離，將計(jì)算結(jié)果與參考解進(jìn)行比較。同時(shí)約定下網(wǎng)格邊界以地下異常體的下邊界為基準(zhǔn)進(jìn)行擴(kuò)展，左右網(wǎng)格邊界以異常體的左右邊界為基準(zhǔn)進(jìn)行擴(kuò)展。約定計(jì)算趨膚深度時(shí)的電阻率為地下電性結(jié)構(gòu)中電阻率最高者。
對(duì)背景介質(zhì)中存在低阻異常體和高阻異常體2種模型進(jìn)行模擬計(jì)算，并在二度體中心上方O點(diǎn)和邊界上方B點(diǎn)同時(shí)進(jìn)行TE模式和TM（transverse magnetic）模式觀測(cè)。
2.2.1 左右網(wǎng)格邊界變化
表3 模型2D－1在O點(diǎn)觀測(cè)結(jié)果Table3 Observed results of model 2D－1on point O
表4 左右網(wǎng)格邊界變化下視電阻率觀測(cè)結(jié)果Table4 Apparent resistivity results when changing left and right mesh boundary
從表4可看出：對(duì)于模型2D－1，TE模式在1δ網(wǎng)格邊界時(shí)，O點(diǎn)和B點(diǎn)觀測(cè)結(jié)果的相對(duì)誤差除250Hz外大部降至1%以下；在2δ的網(wǎng)格邊界時(shí)，O點(diǎn)觀測(cè)結(jié)果的相對(duì)誤差已全降至0.1%以下，B點(diǎn)雖相對(duì)于1δ網(wǎng)格邊界條件時(shí)的相對(duì)誤差有所減小，但部分頻率的誤差仍在0.1%數(shù)量級(jí)上；繼續(xù)增大網(wǎng)格邊界至3δ時(shí)，O點(diǎn)和B點(diǎn)觀測(cè)的相對(duì)誤差已全部降至0.1%以下；繼續(xù)加大網(wǎng)格邊界至5δ和10δ，觀測(cè)結(jié)果與參考結(jié)果相等，誤差為0。
對(duì)于TM模式：在1δ網(wǎng)格邊界時(shí)，計(jì)算的結(jié)果相對(duì)誤差全部在1%以下；增大到2δ網(wǎng)格邊界時(shí)，相對(duì)誤差大部為0；繼續(xù)加大網(wǎng)格邊界，相對(duì)誤差保持為0。對(duì)于模型2D－2，在1δ網(wǎng)格邊界時(shí)，大部分計(jì)算結(jié)果與參考解相等，相對(duì)誤差為0。
值得注意的是，對(duì)模型2D－2進(jìn)行模擬時(shí)，網(wǎng)格邊界變化的量度采用的趨膚深度已相當(dāng)于模型2D－1的3倍。
2.2.2 下網(wǎng)格邊界變化
在對(duì)下邊界研究時(shí)，則把左右網(wǎng)格邊界設(shè)置在“足夠遠(yuǎn)”的距離，以趨膚深度作為量度來改變下網(wǎng)格邊界的邊界距離，計(jì)算結(jié)果并與參考解進(jìn)行比較。仍然以模型2D－1和模型2D－2為計(jì)算模型。模型2D－1的部分計(jì)算結(jié)果如表5。
從表5中看出：根據(jù)TE模式計(jì)算結(jié)果，在1δ網(wǎng)格邊界時(shí)，O點(diǎn)和B點(diǎn)觀測(cè)結(jié)果的相對(duì)誤差已全部降至1%以下；在2δ網(wǎng)格邊界時(shí)，相對(duì)誤差繼續(xù)減小，但大部仍在0.1%的數(shù)量級(jí)上；繼續(xù)增大網(wǎng)格邊界至3δ時(shí)，O點(diǎn)和B點(diǎn)觀測(cè)的相對(duì)誤差已全部降至0.1%以下；加大邊界至5δ和10δ，觀測(cè)結(jié)果與參考結(jié)果相等，誤差為0。
對(duì)于TM模式：在1δ網(wǎng)格邊界時(shí)，計(jì)算的結(jié)果相對(duì)誤差大部在1%以下；增大到2δ網(wǎng)格邊界時(shí)，除B點(diǎn)觀測(cè)的25Hz和250Hz外，大部分結(jié)果的相對(duì)誤差降至0.1%以下；繼續(xù)增大網(wǎng)格邊界，誤差繼續(xù)減少直至為0。
對(duì)于下網(wǎng)格邊界變化下的高阻模型2D－2，計(jì)算的結(jié)果與對(duì)應(yīng)的左右網(wǎng)格邊界的變化結(jié)果相似。即在1δ時(shí)，大部分計(jì)算結(jié)果已與參考解相等，相對(duì)誤差為0。
3 截?cái)噙吔邕x擇
3.1 綜合網(wǎng)格邊界計(jì)算
以上的計(jì)算結(jié)論都是在某一方向的網(wǎng)格邊界位于“足夠遠(yuǎn)處”得到的，現(xiàn)在把以上得出的最佳網(wǎng)格邊界距離綜合在一起進(jìn)行計(jì)算比較。
對(duì)于模型2D－1：計(jì)算TE模式時(shí)，左右網(wǎng)格邊界和下網(wǎng)格邊界都放在3δ處；計(jì)算TM模式時(shí)，則將左右網(wǎng)格邊界和下網(wǎng)格邊界放在2δ處。對(duì)于模型2D－2，把網(wǎng)格邊界全部設(shè)置在1δ處即可。綜合網(wǎng)格邊界條件下部分觀測(cè)結(jié)果如表6所示。
由表6可見，以上利用綜合網(wǎng)格邊界計(jì)算的結(jié)果相對(duì)誤差基本在0.1%以下，說明該網(wǎng)格邊界條件下計(jì)算結(jié)果的可靠性。
3.2 綜合模型計(jì)算
為了進(jìn)一步驗(yàn)證綜合的最佳網(wǎng)格邊界計(jì)算結(jié)果的可靠性，建立1個(gè)較復(fù)雜的二維模型，并采用新的頻率進(jìn)行計(jì)算。模型2D－3（圖4）是在背景介質(zhì)中賦存3種異常體的較復(fù)雜模型，分別在O、O1、O2和B點(diǎn)觀測(cè)，計(jì)算的頻率為0.01、0.10、1.00、10.00、100.00、1 000.00Hz，剖分情況與前面相同。
在模型2D－3中存在多個(gè)異常體，在確定網(wǎng)格邊界變化量度時(shí)，需要幾個(gè)邊界變化量度進(jìn)行比較才能確定。在TE模式計(jì)算時(shí)，對(duì)于低阻體ρ2和ρ4，當(dāng)它們賦存在較高電阻率的背景介質(zhì)中時(shí)，趨膚深度計(jì)算采用背景電阻率，邊界變化量度為3δ＝503；而當(dāng)背景介質(zhì)中存在高阻體ρ3時(shí)，趨膚深度計(jì)算采用高阻體電阻率ρ3，其邊界變化量度為δ ＝503-。為了確保計(jì)算結(jié)果正確和保證精度，網(wǎng)格邊界變化量度應(yīng)為TM模式計(jì)算時(shí)網(wǎng)格邊界變化量度的確定與之類似。表7為模型2D－3的部分計(jì)算結(jié)果。
表5 模型2D－1下網(wǎng)格邊界變化下O點(diǎn)視電阻率觀測(cè)結(jié)果Table5 Apparent resistivity results of model 2D－1on point Owhen changing bottom boundary
表6 綜合網(wǎng)格邊界O點(diǎn)觀測(cè)結(jié)果Table6 Observed results on point Owith compositive boundaries
表7 模型2D－3在O點(diǎn)的觀測(cè)結(jié)果Table7 Observed results of model 2D－3on point O
圖4 模型2D－3Fig.4 Model 2D－3
由表7可見，利用綜合的網(wǎng)格邊界計(jì)算的結(jié)果與參考解的相對(duì)誤差全部在0.1%以下，進(jìn)一步驗(yàn)證了綜合邊界的正確性。
值得注意的是，以上的結(jié)果是在背景電阻率和異常體電阻率最大相差一個(gè)數(shù)量級(jí)的情況下得出的。作者又利用綜合邊界計(jì)算了背景電阻率為100 Ω·m、異常體電阻率分別為1、0.1、0.01、0.001 Ω·m的情況以及背景電阻率為100Ω·m、異常體電阻率為200Ω·m的情況，計(jì)算結(jié)果表明TE模式中誤差最大為0.40%，TM模式誤差基本為0。另外作者還計(jì)算了一些極端的情況，如背景電阻率為1Ω·m、異常體電阻率為2Ω·m的模型和背景電阻率為2Ω·m、異常體電阻率為1Ω·m的模型，發(fā)現(xiàn)參考邊界在此也適用。
4 結(jié)論與建議
利用有限元方法對(duì)大地電磁進(jìn)行模擬計(jì)算時(shí)，由于存在截?cái)噙吔缬绊懀瑥亩绊懹?jì)算結(jié)果和精度，或占用大量內(nèi)存花費(fèi)較長的時(shí)間。筆者通過大量的數(shù)值模擬計(jì)算與對(duì)比，以趨膚深度作為網(wǎng)格邊界變化量度，考慮到大地電磁的精度要求，總結(jié)得出適合二維大地電磁有限元正演模擬計(jì)算的參考網(wǎng)格邊界：
1）利用有限元對(duì)一維地電模型進(jìn)行模擬計(jì)算時(shí)，不存在截?cái)噙吔绲挠绊?，邊界條件自然滿足。
2）對(duì)二維低阻模型進(jìn)行模擬計(jì)算時(shí)，以背景電阻率計(jì)算趨膚深度作為網(wǎng)格邊界變化量度。對(duì)于TE模式，當(dāng)左右網(wǎng)格邊界和下網(wǎng)格邊界都放在3δ處時(shí)，截?cái)噙吔绲挠绊懸芽珊雎浴?duì)于TM模式，左右網(wǎng)格邊界和下網(wǎng)格邊界在2δ處，截?cái)噙吔缬绊懣梢院雎浴?/p>
3）對(duì)二維高阻模型進(jìn)行模擬計(jì)算時(shí)，以高阻體電阻率計(jì)算趨膚深度作為網(wǎng)格邊界變化量度。不管TE模式還是TM模式，當(dāng)左右網(wǎng)格邊界和下網(wǎng)格邊界位于1δ時(shí)，截?cái)噙吔缬绊懸芽珊雎浴?/p>
4）當(dāng)二維地電模型中，同時(shí)存在低阻異常體和高阻異常體時(shí)，分別以背景介質(zhì)電阻率和高阻體中電阻率最高者計(jì)算趨膚深度，然后根據(jù)參考網(wǎng)格邊界選擇最佳的網(wǎng)格邊界。
（References）：
［1］Coggon J H.Electromagnetic and Electrical Modeling by Finite Element Method［J］.Geophysics，1971，36（1）：132－155.
［2］Rodi W L.A Technique for Improving the Accuracy of Finite Element Solution for Magnetotelluric Data［J］.Geophysics，1976，44（2）：483－506.
［3］Wannamaker P E，Stodt J A，Rijo W L.Two－Dimensional Topographic Responses in Magnetotelluric Modeling Using Finite Element［J］.Geophysics，1986，51（11）：2131－2144.
［4］Wannamaker P E，Stodt J A，Rijo W L.A Stable Finite Element Solution for Two－Dimensional Magnetotelluric Modelling［J］.Geophysics，1987，88（1）：277－296.
［5］陳樂壽.有限元法在大地電磁場(chǎng)正演計(jì)算中的應(yīng)用與改進(jìn)［J］.石油物探，1981，20（3）：84－103.Chen Leshou.The Improvement and Application of Finite Element Method to 2DForward Solution of Magnetotelluric Sounding［J］.Geophysical Prospecting for Petroleum，1981，20（3）：84－103.
［6］胡建德，蔡剛.用三角形二次插值有限元法計(jì)算二維大地電磁測(cè)深曲線［J］.石油地球物理勘探，1984，19（4）：358－367.Hu Jiande，Cai Gang.Calculating 2DMagnetotelluric Sounding Curve with the Twice Interpolation Triangle Finite Element Method ［J］. Oil Geophysical Prospecting，1984，19（4）：358－367.
［7］史明娟，徐世浙，劉斌.大地電磁二次函數(shù)插值的有限元法正演模擬［J］.地球物理學(xué)報(bào)，1997，40（3）：421－430.Shi Mingjuan，Xu Shizhe，Liu Bin.Finite Element Method Using Quadratic Element in MT Forward Modeling［J］.Acta Geophysica Sinica，1997，40（3）：421－430.
［8］楊長福.大地電磁二維對(duì)稱各向異性介質(zhì)的有限元數(shù)值模擬［J］.西北地震學(xué)報(bào)，1997，19（2）：27－33.Yang Changfu.MT Numerical Simulation of Symmetrically 2DAnisotropic Media Based on the Finite Element Method ［J］. Northwestern Seismological Journal，1997，19（2）：27－33.
［9］陳小斌.有限元直接迭代算法［J］.物探化探計(jì)算技術(shù)，1999，21（2）：165－171.Chen Xiaobin.Direct Iterative Algorithm of Finite Element［J］.Computing Techniques for Geophysical and Geochemical Exploration，1999，21（2）：165－171.
［10］陳樂壽，劉任，王天生.大地電磁測(cè)深資料處理與解釋［M］.北京：石油工業(yè)出版社，1989：1－160.Chen Leshou， Liu Ren， Wang Tiansheng.Magnetotelluric Sounding Data Processing and Interpretation［M］.Beijing：Petroleum Industry Press，1989：1－160.
［11］徐世浙.地球物理中的有限單元法［M］.北京：科學(xué)出版社，1994：220－241.Xu Shizhe.The Finite Element in Geophysics［M］.Beijing：Science Press，1994：220－241.
［12］趙廣茂，李桐林，王大勇，等.基于二次場(chǎng)二維起伏地形MT有限元數(shù)值模擬［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2008，38（6）：1055－1059.Zhao Guangmao，Li Tonglin，Wang Dayong，et al.Secondary Field－Based Two－Dimensional Topographic Numerical Simulation in Magnetoteilurics by Finite Element Method［J］.Journal of Jilin University：Earth Science Edition，2008，38（6）：1055－1059.
［13］劉長生，湯井田，任政勇，等.基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬［J］.中南大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，41（5）：1855－1859.Liu Changsheng，Tang Jingtian，Ren Zhengyong，et al.Three－Dimension Magnetotellurics Modeling by Adaptive Edge Finite－Element Using Unstructured Meshes［J］.Journal of Central South University：Science and Technology，2010，41（5）：1855－1859.
［14］Sharma P S，Kaikkonen P.An Automated Finite Element Mesh Generation and Element Coding in 2－D Electromagnetic Inversion［J］.Geophysics，1998，34（3）：93－114.
猜你喜歡
應(yīng)用數(shù)學(xué)(2020年4期)2020-12-28 00:36:38
數(shù)學(xué)物理學(xué)報(bào)(2020年5期)2020-11-26 06:06:30
甘肅科技(2020年20期)2020-04-13 00:30:40
水利科技與經(jīng)濟(jì)(2016年3期)2016-04-22 01:05:00
西南石油大學(xué)學(xué)報(bào)(自然科學(xué)版)(2015年4期)2015-08-20 09:05:16
物探化探計(jì)算技術(shù)(2015年2期)2015-02-28 17:42:42
四川師范大學(xué)學(xué)報(bào)（自然科學(xué)版）(2015年1期)2015-02-28 14:07:21
火炸藥學(xué)報(bào)(2014年3期)2014-03-20 13:17:39
建筑材料學(xué)報(bào)(2014年3期)2014-03-11 17:08:01
東北師大學(xué)報(bào)（自然科學(xué)版）(2014年1期)2014-02-27 08:02:09
