一種各向異性介質(zhì)大地電磁無(wú)網(wǎng)格數(shù)值模擬方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及地球物理勘探領(lǐng)域中地下介質(zhì)的各向異性問(wèn)題，提出了一種基于離散 節(jié)點(diǎn)構(gòu)造形函數(shù)的無(wú)網(wǎng)格數(shù)值模擬方法，適用于大地電磁等地球物理勘探方法。
【背景技術(shù)】
[0002] 地球內(nèi)部的結(jié)構(gòu)和屬性是地球物理學(xué)研究的核心內(nèi)容。21世紀(jì)的地球科學(xué)是介 質(zhì)，構(gòu)造和深層動(dòng)力學(xué)過(guò)程的橫向不均勻性和各向異性的時(shí)代。隨著現(xiàn)代觀測(cè)技術(shù)的不斷 進(jìn)步以及認(rèn)識(shí)水平的不斷提高，各向異性問(wèn)題逐漸引起了人們的廣泛重視，成為地球物理 學(xué)研究的熱點(diǎn)。
[0003] 大地電磁測(cè)深法(Magnetotelluric Sounding，MT)是利用天然交變電磁場(chǎng)研究地 球電性結(jié)構(gòu)的一種地球物理勘探方法。該方法具有對(duì)地下高導(dǎo)構(gòu)造反映敏感的特點(diǎn)，被廣 泛應(yīng)用于地殼、上地幔的研究中?，F(xiàn)有的大地電磁正演計(jì)算方法主要有積分方程法（IEM)， 有限差分法(FDM)，有限單元法(FEM)等。
[0004] CN201410370018.0公開(kāi)了一種時(shí)移大地電磁信號(hào)采集和反演方法，該方法消除了 時(shí)移數(shù)據(jù)中的系統(tǒng)誤差和正演過(guò)程中的計(jì)算誤差，減小了系統(tǒng)誤差對(duì)于所得的不同時(shí)刻模 型差異的影響，使得反演所得的模型變化更接近真實(shí)情況。
[0005] CN201010597160.0公開(kāi)了一種各向異性三維疊前時(shí)間偏移方法，該方法考慮地球 介質(zhì)速度各向異性對(duì)地震波傳播的走時(shí)和幅值的影響，能在偏移過(guò)程中自主決定偏移速度 和各向異性參數(shù)，該方法主要用在于地震勘探中反射地震資料處理中。
[0006] 以上方法都是基于網(wǎng)格實(shí)現(xiàn)的，在計(jì)算實(shí)現(xiàn)的過(guò)程中存在著生成網(wǎng)格成本高，場(chǎng) 值變化劇烈的地方精度低，自適應(yīng)分析困難等缺點(diǎn)。本發(fā)明針對(duì)以上的不足，提出用無(wú)網(wǎng)格 徑向基點(diǎn)插值法(Radical Point Interpolation Method，RPIM)模擬大地電磁的各向異性 問(wèn)題，可以實(shí)現(xiàn)復(fù)雜模型的正演計(jì)算。

【發(fā)明內(nèi)容】

[0007] 本發(fā)明所要解決的問(wèn)題在于提供一種針對(duì)地下空間各向異性介質(zhì)的，不依賴于網(wǎng) 格構(gòu)造形函數(shù)，能夠進(jìn)行復(fù)雜介質(zhì)大地電磁正演數(shù)值模擬方法。
[0008] 本發(fā)明是這樣實(shí)現(xiàn)的，一種各向異性介質(zhì)大地電磁無(wú)網(wǎng)格數(shù)值模擬方法，包括如 下的步驟：
[0009] 1)從各向異性大地電磁邊值問(wèn)題出發(fā)，構(gòu)造等價(jià)泛函，推導(dǎo)對(duì)應(yīng)無(wú)網(wǎng)格徑向基點(diǎn) 插值法的等價(jià)線性方程組；
[0010] 2)讀取當(dāng)前模型參數(shù)，包括頻率參數(shù)，節(jié)點(diǎn)坐標(biāo)，背景單元，支持域，形狀參數(shù)，極 化模式等；
[0011] 3)對(duì)當(dāng)前極化模式進(jìn)行判斷，若為TE極化模式，則計(jì)算區(qū)域包含空氣層，若為TM極 化模式，則不含空氣層；
[0012] 4)對(duì)所有背景網(wǎng)格進(jìn)行循環(huán)，對(duì)背景網(wǎng)格的所有高斯積分點(diǎn)循環(huán)，搜索該高斯積 分點(diǎn)支持域內(nèi)的有效節(jié)點(diǎn)，計(jì)算支持域內(nèi)節(jié)點(diǎn)處的形函數(shù)，求取系數(shù)矩陣和右端項(xiàng)；
[0013] 5)加載本質(zhì)邊界條件，求解線性方程組，得到各個(gè)節(jié)點(diǎn)的場(chǎng)值；
[0014] 6)由視電阻率計(jì)算公式，代入場(chǎng)值求取地面處各個(gè)方向的視電阻率和相位。
[0015] 進(jìn)一步地，步驟1中，各向異性大地電磁邊值問(wèn)題，假設(shè)一個(gè)電性主軸垂直于層面 方面，另外一個(gè)主軸平行于層面方向，構(gòu)造各向異性介質(zhì)模型的電導(dǎo)率張量模型。
[0016] 進(jìn)一步地，步驟4中，計(jì)算支持域內(nèi)節(jié)點(diǎn)處的形函數(shù)，形函數(shù)采用復(fù)合2次徑向基函 數(shù)(MQ-RBF)構(gòu)造。
[0017]進(jìn)一步地，步驟5中，線性方程組的求解采用Krylov子空間的正則化擬殘量極小方 法(Quasi-minimal Residual method，QMR)，實(shí)現(xiàn)了大型稀疏線性方程組高效，精確的求 解。
[0018] 本發(fā)明與現(xiàn)有技術(shù)相比，有益效果在于:本發(fā)明針對(duì)實(shí)際大地電磁探測(cè)地下介質(zhì) 中廣泛存在的各向異性問(wèn)題，假設(shè)一個(gè)電性主軸垂直于層面方面，另外一個(gè)主軸平行于層 面方向，構(gòu)造了各向異性介質(zhì)模型的電導(dǎo)率張量模型。針對(duì)傳統(tǒng)的正演模擬方法依賴于網(wǎng) 格，物性參數(shù)復(fù)雜分布適應(yīng)性差的缺點(diǎn)，提出用無(wú)網(wǎng)格徑向基點(diǎn)插值法模擬大地電磁的各 向異性問(wèn)題，可以有效地克服單純多項(xiàng)式基點(diǎn)插值方法中存在的奇異性問(wèn)題，形函數(shù)光滑 穩(wěn)定，實(shí)現(xiàn)了電磁法高精度，自適應(yīng)的數(shù)值模擬。
【附圖說(shuō)明】
[0019] 圖1是本發(fā)明實(shí)施例提供的無(wú)網(wǎng)格法場(chǎng)節(jié)點(diǎn)、高斯點(diǎn)、求解域及其邊界，支持域與 背景網(wǎng)格示意圖；
[0020] 圖2是本發(fā)明實(shí)施例提供的各向異性介質(zhì)的電性主軸與地下空間坐標(biāo)系示意圖；
[0021] 圖3是本發(fā)明實(shí)施例提供的各向異性介質(zhì)大地電磁無(wú)網(wǎng)格徑向基點(diǎn)插值法數(shù)值模 擬流程圖；
[0022]圖4是本發(fā)明實(shí)施例建立的層狀各向異性介質(zhì)模型，其中旋轉(zhuǎn)角度為30° ;
[0023] 圖5是本發(fā)明實(shí)施例各項(xiàng)異性層狀模型擬解析解與無(wú)網(wǎng)格解的對(duì)比；
[0024] 圖6是本發(fā)明實(shí)施例建立的均勻介質(zhì)中含有各向異性異常體模型；
[0025] 圖7是本發(fā)明實(shí)施例建立的均勻介質(zhì)中含有各向異性異常體模型各個(gè)方向視電阻 率計(jì)算結(jié)果。
【具體實(shí)施方式】
[0026] 為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點(diǎn)更加清楚明白，以下結(jié)合實(shí)施例，對(duì)本發(fā)明 進(jìn)行進(jìn)一步詳細(xì)說(shuō)明。應(yīng)當(dāng)理解，此處所描述的具體實(shí)施例僅僅用以解釋本發(fā)明，并不用于 限定本發(fā)明。
[0027] 參見(jiàn)圖1，本發(fā)明實(shí)施例提供的無(wú)網(wǎng)格法計(jì)算點(diǎn)1、求解域2,背景單元3,場(chǎng)節(jié)點(diǎn)4， 支持域5與求解域邊界6示意圖。
[0028] 求解域是整個(gè)邊值問(wèn)題所要求解的區(qū)域，支持域是構(gòu)造無(wú)網(wǎng)格法形函數(shù)所選擇的 區(qū)域，計(jì)算點(diǎn)是支持域的中心。在計(jì)算中，計(jì)算點(diǎn)常常選擇為高斯積分點(diǎn)。背景單元是便于 進(jìn)行域積分而劃分的區(qū)域，背景單元并不依賴節(jié)點(diǎn)而存在。
[0029] 參見(jiàn)圖3,各向異性介質(zhì)大地電磁無(wú)網(wǎng)格數(shù)值模擬方法流程圖，包括如下的步驟：
[0030] 1)從各向異性大地電磁邊值問(wèn)題出發(fā)，構(gòu)造等價(jià)泛函，推導(dǎo)對(duì)應(yīng)無(wú)網(wǎng)格徑向基點(diǎn) 插值法的等價(jià)線性方程組；
[0031 ] 2)讀取當(dāng)前模型參數(shù)，包括頻率參數(shù)，節(jié)點(diǎn)坐標(biāo)，背景單元，支持域，形狀參數(shù)，極 化模式等；
[0032] 3)對(duì)當(dāng)前極化模式進(jìn)行判斷，若為TE極化模式，則計(jì)算區(qū)域包含空氣層，若為TM極 化模式，則不含空氣層；
[0033] 4)對(duì)所有背景網(wǎng)格進(jìn)行循環(huán)，對(duì)背景網(wǎng)格的所有高斯積分點(diǎn)循環(huán)，搜索該高斯積 分點(diǎn)支持域內(nèi)的有效節(jié)點(diǎn)，計(jì)算支持域內(nèi)節(jié)點(diǎn)處的形函數(shù)，求取系數(shù)矩陣和右端項(xiàng)；
[0034] 5)加載本質(zhì)邊界條件，求解線性方程組，得到各個(gè)節(jié)點(diǎn)的場(chǎng)值；
[0035] 6)由視電阻率計(jì)算公式，代入場(chǎng)值求取地面處各個(gè)方向的視電阻率和相位。
[0036] 當(dāng)步驟4中計(jì)算系數(shù)矩陣和右端項(xiàng)結(jié)束后檢測(cè)高斯積分點(diǎn)循環(huán)是否結(jié)束，結(jié)束后 進(jìn)一步檢測(cè)背景網(wǎng)格是否循環(huán)完畢，高斯積分點(diǎn)循環(huán)未結(jié)束時(shí)，返回步驟4對(duì)高斯積分點(diǎn)循 環(huán)，背景網(wǎng)格循環(huán)完畢后進(jìn)入下一步，當(dāng)背景網(wǎng)格未循環(huán)完畢時(shí)，則返回到步驟4對(duì)背景網(wǎng) 絡(luò)進(jìn)行循環(huán)。
[0037]步驟1中各向異性大地電磁邊值問(wèn)題 [0038]在介質(zhì)中，頻率域電磁場(chǎng)方程是： VxE -iconII /
[0039] '' 、 ( 1 ) V x // = ^a-ico〇) E
[0040] 在地下介質(zhì)空間中，對(duì)大地電磁探測(cè)所涉及的頻率而言，電導(dǎo)率〇>> e，因而也可 不考慮e的各向異性問(wèn)題。電導(dǎo)率〇在各向同性介質(zhì)中是標(biāo)量，在各向異性介質(zhì)中是張量。本 發(fā)明只考慮〇的各向異性問(wèn)題。
[0041] 如圖2所示，建立如圖所示的坐標(biāo)系，設(shè)平行層面的電導(dǎo)率為〇//，垂直層面的電導(dǎo) 率為〇丄，在層面坐標(biāo)系x ' y ' z '中x '平行層面，y '垂直層面，z '平行走向，電導(dǎo)率張量為
(2)
[0043] 在地面坐標(biāo)系中，取走向?yàn)閦軸（與z'平行），x軸與z軸垂直，保持水平，y軸垂直向 上。在坐標(biāo)系xyz中，電導(dǎo)率〇張量為
[0044] o=A〇，At (3)
[0045] 其中A為坐標(biāo)變換張量
(4)
[0047]其中a是x'軸與地下空間坐標(biāo)系x軸的夾角。
[0048]將（1)式按分量展開(kāi)，并考慮到則將電導(dǎo)率張量代入，引入二維算子
經(jīng)化簡(jiǎn)可得各向異性大地電磁邊值問(wèn)題如下： (5)
[0050]其中U指的是場(chǎng)函數(shù)Ez或者Hz，Q指的是求解域，AB是求解域的上邊界，AC和BD指的 是求解域左右邊界，CD指的是求解域下邊界，式中：
《代表圓頻率，y代表介 質(zhì)的磁導(dǎo)率，在TE極化模式下，u = Ez，A = 〇//-i ? e :。在TM極化模式下，u = Hz，A=i ? y，T的表達(dá)式如下：

[0053 ]步驟4中支持域內(nèi)徑向基點(diǎn)插值法建立形函數(shù)
[0054] 如圖1所示，在求解域Q內(nèi)，徑向基點(diǎn)插值函數(shù)(RPIM)可表示為
[0056]式中Ri(X)為徑向基函數(shù)(RBF)，n為RBFs的個(gè)數(shù)，Pj(X)為空間坐標(biāo)XT=(x，y)中的 單項(xiàng)式，m為多項(xiàng)式基函數(shù)的個(gè)數(shù)。&1和匕為待定常數(shù)。本發(fā)明選用的徑向基函數(shù)為復(fù)合2次 徑向基函數(shù)(MQ-RBF)，其表達(dá)式為：
[0057] Ri(x,y) = (ri2+(acdc)2)q(ac 0) (9)
[0058] 為確定式(8)中的ai和bj，需形成計(jì)算點(diǎn)X的支持域，其中包括n個(gè)場(chǎng)節(jié)點(diǎn)。使式(8) 滿足計(jì)算點(diǎn)X周圍的n個(gè)節(jié)點(diǎn)值以確定系數(shù) &1和匕，這將