一種隨機(jī)導(dǎo)電媒質(zhì)模型的大地電磁無(wú)網(wǎng)格數(shù)值模擬方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及地球物理勘探領(lǐng)域中大地電磁法的數(shù)值模擬方法，在計(jì)算過(guò)程中脫離 了傳統(tǒng)方法中對(duì)于網(wǎng)格的依賴(lài)，適用于隨機(jī)導(dǎo)電媒質(zhì)模型等復(fù)雜模型。
【背景技術(shù)】
[0002] 大地電磁測(cè)深法（Magnetotelluric Sounding)，簡(jiǎn)稱(chēng)MT，是蘇聯(lián)學(xué)者 Tikkonov(1950)和法國(guó)學(xué)者Cagniard(1953)50年代初提出來(lái)的、利用天然交變電磁場(chǎng)研 究地球電性結(jié)構(gòu)的一種地球物理勘探方法。由于它不用人工供電，工作方便，不受高阻層屏 蔽，勘探深度大，勘探成本低，因此受到廣泛的應(yīng)用。
[0003] 在建立模型進(jìn)行研究時(shí)，通常將地下的介質(zhì)看作均勻的或者是層狀介質(zhì)，但在實(shí) 際地球物理勘探中，地層具有復(fù)雜的非均質(zhì)特征，由于在介質(zhì)中的這些干擾對(duì)整體影響太 小而常常被忽略，甚至將這種小尺度上的非均勻性所形成的各種復(fù)雜的干擾作為噪聲而處 理掉，這樣不利于準(zhǔn)確地描述真實(shí)的地下空間構(gòu)造特征。使用隨機(jī)導(dǎo)電媒質(zhì)模型能夠更加 準(zhǔn)確地描述地下媒質(zhì)的非均質(zhì)性，建立的模型能夠更加真實(shí)地反映地層的情況，使研究具 有更好的實(shí)用價(jià)值。
[0004] 數(shù)值模擬的正演計(jì)算是地球物理反演的基礎(chǔ)，對(duì)提高地球物理勘探的效果具有重 要的意義。現(xiàn)有的大地電磁正演計(jì)算方法主要有積分方程法（IEM)，有限差分法（FDM)，有 限單元法（FEM)等。這些方法都是基于網(wǎng)格實(shí)現(xiàn)的，在計(jì)算實(shí)現(xiàn)的過(guò)程中存在著生成網(wǎng)格 成本高，場(chǎng)值變化劇烈的地方精度低，自適應(yīng)分析困難等缺點(diǎn)。

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

[0005] 本發(fā)明所要解決的技術(shù)問(wèn)題在于提供一種針對(duì)地下空間非均質(zhì)特征的，基于離散 節(jié)點(diǎn)構(gòu)造形函數(shù)，擺脫對(duì)于網(wǎng)格依賴(lài)的復(fù)雜介質(zhì)正演數(shù)值模擬方法。
[0006] 本發(fā)明是這樣實(shí)現(xiàn)的，一種隨機(jī)導(dǎo)電媒質(zhì)模型的大地電磁無(wú)網(wǎng)格數(shù)值模擬方法， 包括如下的步驟：
[0007] 1)構(gòu)造隨機(jī)導(dǎo)電媒質(zhì)模型的電性參數(shù)；
[0008] 2)設(shè)置模型參數(shù)，包括頻率參數(shù)，節(jié)點(diǎn)坐標(biāo)，背景網(wǎng)格，支持域等；
[0009] 3)對(duì)單個(gè)頻率進(jìn)行循環(huán)，對(duì)所有背景網(wǎng)格進(jìn)行循環(huán)，包括：
[0010] 對(duì)背景網(wǎng)格的所有高斯積分點(diǎn)循環(huán)，搜索該高斯積分點(diǎn)支持域的有效節(jié)點(diǎn)；
[0011] 計(jì)算支持域內(nèi)節(jié)點(diǎn)處的形函數(shù)；
[0012] 4)建立和邊值問(wèn)題的等價(jià)泛函，加載本質(zhì)邊界條件、混合邊界條件，分別得到系數(shù) 矩陣和右端項(xiàng)，組裝到線(xiàn)性方程組中；
[0013] 5)求解線(xiàn)性方程組，得到各個(gè)節(jié)點(diǎn)的場(chǎng)值；
[0014] 6)由卡尼亞視電阻率計(jì)算公式，代入場(chǎng)值求取地面處的視電阻率和相位。
[0015] 進(jìn)一步地，步驟1中，隨機(jī)導(dǎo)電媒質(zhì)模型的電性參數(shù)采用均值為零的橢圓自相關(guān) 函數(shù)的二階平穩(wěn)隨機(jī)過(guò)程來(lái)構(gòu)造。
[0016] 進(jìn)一步地，步驟4中，本質(zhì)邊界條件采用拉格朗日乘子法加載，構(gòu)造了引入約束條 件的修正泛函，實(shí)現(xiàn)了本質(zhì)邊界條件的精確處理。
[0017] 進(jìn)一步地，步驟5中，線(xiàn)性方程組的求解采用基于雙邊Lanczos算法和殘差正交子 空間的穩(wěn)定雙共輒梯度（BICGSTAB)迭代算法，采用不完全LU分解方法預(yù)處理系數(shù)矩陣，實(shí) 現(xiàn)了線(xiàn)性方程快速，精確的求解。
[0018] 本發(fā)明與現(xiàn)有技術(shù)相比，有益效果在于：本發(fā)明針對(duì)實(shí)際大地電磁探測(cè)中地下介 質(zhì)中存在的大量分布不規(guī)則的微小異常造成的不相干擾動(dòng)的情況，基于隨機(jī)建模過(guò)程的譜 分解理論和混合型自相關(guān)函數(shù)理論構(gòu)造地下空間的電性參數(shù)模型。針對(duì)傳統(tǒng)的網(wǎng)格化正演 方式的不足，提出了無(wú)網(wǎng)格法模擬的思想及方案。由于無(wú)網(wǎng)格法只需要節(jié)點(diǎn)離散場(chǎng)值，克服 了網(wǎng)格化正演方法的不足，介質(zhì)物性參數(shù)的加載非常方便，非常適合于隨機(jī)導(dǎo)電媒質(zhì)等復(fù) 雜模型的計(jì)算。通過(guò)對(duì)比表明，隨機(jī)導(dǎo)電媒質(zhì)模型能夠更加真實(shí)有效地反應(yīng)地下實(shí)際介質(zhì) 的分布情況，正演計(jì)算剖面與實(shí)際剖面情況更加相符，從而更加有利于認(rèn)識(shí)大地電磁波的 傳播規(guī)律，指導(dǎo)大地電磁法的數(shù)據(jù)解釋。
【附圖說(shuō)明】
[0019] 圖1是本發(fā)明實(shí)施例提供的無(wú)網(wǎng)格法計(jì)算點(diǎn)、場(chǎng)節(jié)點(diǎn)、求解域，支持域與背景網(wǎng)格 示意圖；
[0020] 圖2是本發(fā)明實(shí)施例提供的在二維空間中大地電磁邊值問(wèn)題示意圖；
[0021] 圖3是本發(fā)明實(shí)施例提供的隨機(jī)導(dǎo)電媒質(zhì)模型大地電磁無(wú)網(wǎng)格發(fā)正演模擬方法 流程圖；
[0022] 圖4是本發(fā)明實(shí)施例建立的層狀介質(zhì)模型；
[0023] 圖5是本發(fā)明實(shí)施例層狀模型有限元解，解析解與無(wú)網(wǎng)格解的對(duì)比；
[0024] 圖6是本發(fā)明實(shí)施例建立的二維大地隨機(jī)導(dǎo)電媒質(zhì)模型；
[0025] 圖7是本發(fā)明實(shí)施例無(wú)網(wǎng)格法求解大地電磁二維大地隨機(jī)導(dǎo)電媒質(zhì)模型視電阻 率結(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)4、場(chǎng)節(jié)點(diǎn)2、求解域1，支持域3與 背景網(wǎng)格5示意圖。
[0028] 求解域是整個(gè)邊值問(wèn)題所要求解的區(qū)域，支持域是構(gòu)造無(wú)網(wǎng)格法形函數(shù)所選擇的 區(qū)域，計(jì)算點(diǎn)是支持域的中心。在計(jì)算中，計(jì)算點(diǎn)常常選擇為高斯積分點(diǎn)。背景網(wǎng)格是進(jìn)行 域積分而劃分的區(qū)域。
[0029] 參見(jiàn)圖3,隨機(jī)導(dǎo)電媒質(zhì)模型大地電磁無(wú)網(wǎng)格正演模擬方法流程圖，包括如下的步 驟：
[0030] 1)構(gòu)造隨機(jī)導(dǎo)電媒質(zhì)模型的電性參數(shù)；
[0031] 2)設(shè)置模型參數(shù)，包括頻率參數(shù)，節(jié)點(diǎn)坐標(biāo)，背景網(wǎng)格，支持域等；
[0032] 3)對(duì)單個(gè)頻率進(jìn)行循環(huán)，對(duì)所有背景網(wǎng)格進(jìn)行循環(huán)，包括：
[0033] 對(duì)背景網(wǎng)格的所有高斯積分點(diǎn)循環(huán)，搜索該高斯積分點(diǎn)支持域的有效節(jié)點(diǎn)；
[0034] 計(jì)算支持域內(nèi)節(jié)點(diǎn)處的形函數(shù)；
[0035] 4)利用拉格朗日乘子法處理上邊界，建立和邊值問(wèn)題的等價(jià)泛函，加載本質(zhì)邊界 條件、混合邊界條件，分別得到系數(shù)矩陣和右端項(xiàng)，組裝到線(xiàn)性方程組中；
[0036] 5)采用不完全LU分解預(yù)處理的穩(wěn)定雙共輒梯度法求解線(xiàn)性方程組，得到各個(gè)節(jié) 點(diǎn)的場(chǎng)值；
[0037] 6)由卡尼亞視電阻率計(jì)算公式，代入場(chǎng)值求取地面處的視電阻率和相位。
[0038] 當(dāng)步驟3中計(jì)算系數(shù)矩陣結(jié)束后檢測(cè)高斯積分點(diǎn)循環(huán)是否結(jié)束，結(jié)束后進(jìn)行下一 步檢測(cè)背景網(wǎng)格是否循環(huán)完畢，高斯積分點(diǎn)循環(huán)未結(jié)束時(shí)，返回步驟3對(duì)高斯積分點(diǎn)循環(huán)， 背景網(wǎng)格循環(huán)完畢后進(jìn)入下一步，當(dāng)背景網(wǎng)格未循環(huán)完畢，則返回到步驟3對(duì)背景網(wǎng)絡(luò)進(jìn) 行循環(huán)。
[0039] 步驟1中構(gòu)造二階平穩(wěn)過(guò)程m(x)的隨機(jī)導(dǎo)電媒質(zhì)模型：
[0040] 隨機(jī)導(dǎo)電媒質(zhì)模型主要由非均勻性大、小兩種尺度組成。大尺度描述介質(zhì)的背景 特性，而小尺度則描述加在背景模型上的隨機(jī)擾動(dòng)。隨機(jī)導(dǎo)電媒質(zhì)模型通常用一個(gè)均值為 零的二階平穩(wěn)隨機(jī)過(guò)程來(lái)描述。在非劇烈變化的地質(zhì)條件下，用m(x)表示非均勻介質(zhì)，mi(x)表示均勻介質(zhì)（介質(zhì)背景），Sm(x)表示小尺度的隨機(jī)擾動(dòng)。采用如下公式描述隨機(jī) 導(dǎo)電媒質(zhì)模型：
[0042] (1)式中，f(x)表示介質(zhì)的相對(duì)擾動(dòng)特性，且
>x表示空間位置矢量。 通?？梢杂脵E圓自相關(guān)函數(shù)來(lái)構(gòu)造隨機(jī)導(dǎo)電媒質(zhì)模型，其表達(dá)式如下：
[0044] ⑵式中，r表示模糊度因子（當(dāng)r = 0時(shí)，為高斯型橢圓自相關(guān)函數(shù)，當(dāng)r = 1時(shí)， 為指數(shù)型橢圓自相關(guān)函數(shù)，r介于兩者之間時(shí)，則為混合型自相關(guān)函數(shù)），a、b分別表示X、z 方向的自相關(guān)長(zhǎng)度，對(duì)式（2)做二維傅里葉變換，得到F(k)，再利用F(k)構(gòu)造能量譜密度 函數(shù)：
[0046] (3)式中，k = (kx，kz)T，用隨機(jī)數(shù)發(fā)生器在區(qū)間[0· 2 π ]上生成服從均勻分布且獨(dú) 立的二維隨機(jī)相位切斤），對(duì)能量譜密度函數(shù)R(k)加入隨機(jī)相位經(jīng)變換后得到：
[0048] 式中，W(k)為窗口函數(shù)，其目的是減小自相關(guān)函數(shù)的低頻能量，對(duì)（4)式進(jìn)行二維 傅里葉逆變換，得到f'（X，Z)，計(jì)算f'（X，Z)的均值和方差即：
[0051] 再將式（5)、（6)進(jìn)行規(guī)范化到均值為零、方差為ε2,并以f(x，z)為自相關(guān)函數(shù)的 二階平穩(wěn)隨機(jī)過(guò)程：
[0053] 通過(guò)上述計(jì)算就可以得到地下空間介質(zhì)的隨機(jī)擾動(dòng)，代入（1)式構(gòu)造相應(yīng)尺度的 二維隨機(jī)導(dǎo)電媒質(zhì)模型。隨機(jī)導(dǎo)電媒質(zhì)的介電常數(shù)和電導(dǎo)率均按照上述方法構(gòu)建。
[0054] 步驟3中支持域內(nèi)滑動(dòng)最小二乘法建立形函數(shù)：
[0055] 如圖1所示，在求解域Ω內(nèi)，假設(shè)場(chǎng)函數(shù)u(x)在N個(gè)節(jié)點(diǎn)上的場(chǎng)值u已知，即：