一種各向異性介質大地電磁無網格數值模擬方法
【技術領域】
[0001] 本發明涉及地球物理勘探領域中地下介質的各向異性問題，提出了一種基于離散 節點構造形函數的無網格數值模擬方法，適用于大地電磁等地球物理勘探方法。
【背景技術】
[0002] 地球內部的結構和屬性是地球物理學研究的核心內容。21世紀的地球科學是介 質，構造和深層動力學過程的橫向不均勻性和各向異性的時代。隨著現代觀測技術的不斷 進步以及認識水平的不斷提高，各向異性問題逐漸引起了人們的廣泛重視，成為地球物理 學研究的熱點。
[0003] 大地電磁測深法(Magnetotelluric Sounding，MT)是利用天然交變電磁場研究地 球電性結構的一種地球物理勘探方法。該方法具有對地下高導構造反映敏感的特點，被廣 泛應用于地殼、上地幔的研究中。現有的大地電磁正演計算方法主要有積分方程法（IEM)， 有限差分法(FDM)，有限單元法(FEM)等。
[0004] CN201410370018.0公開了一種時移大地電磁信號采集和反演方法，該方法消除了 時移數據中的系統誤差和正演過程中的計算誤差，減小了系統誤差對于所得的不同時刻模 型差異的影響，使得反演所得的模型變化更接近真實情況。
[0005] CN201010597160.0公開了一種各向異性三維疊前時間偏移方法，該方法考慮地球 介質速度各向異性對地震波傳播的走時和幅值的影響，能在偏移過程中自主決定偏移速度 和各向異性參數，該方法主要用在于地震勘探中反射地震資料處理中。
[0006] 以上方法都是基于網格實現的，在計算實現的過程中存在著生成網格成本高，場 值變化劇烈的地方精度低，自適應分析困難等缺點。本發明針對以上的不足，提出用無網格 徑向基點插值法(Radical Point Interpolation Method，RPIM)模擬大地電磁的各向異性 問題，可以實現復雜模型的正演計算。
【發明內容】
[0007] 本發明所要解決的問題在于提供一種針對地下空間各向異性介質的，不依賴于網 格構造形函數，能夠進行復雜介質大地電磁正演數值模擬方法。
[0008] 本發明是這樣實現的，一種各向異性介質大地電磁無網格數值模擬方法，包括如 下的步驟：
[0009] 1)從各向異性大地電磁邊值問題出發，構造等價泛函，推導對應無網格徑向基點 插值法的等價線性方程組；
[0010] 2)讀取當前模型參數，包括頻率參數，節點坐標，背景單元，支持域，形狀參數，極 化模式等；
[0011] 3)對當前極化模式進行判斷，若為TE極化模式，則計算區域包含空氣層，若為TM極 化模式，則不含空氣層；
[0012] 4)對所有背景網格進行循環，對背景網格的所有高斯積分點循環，搜索該高斯積 分點支持域內的有效節點，計算支持域內節點處的形函數，求取系數矩陣和右端項；
[0013] 5)加載本質邊界條件，求解線性方程組，得到各個節點的場值；
[0014] 6)由視電阻率計算公式，代入場值求取地面處各個方向的視電阻率和相位。
[0015] 進一步地，步驟1中，各向異性大地電磁邊值問題，假設一個電性主軸垂直于層面 方面，另外一個主軸平行于層面方向，構造各向異性介質模型的電導率張量模型。
[0016] 進一步地，步驟4中，計算支持域內節點處的形函數，形函數采用復合2次徑向基函 數(MQ-RBF)構造。
[0017]進一步地，步驟5中，線性方程組的求解采用Krylov子空間的正則化擬殘量極小方 法(Quasi-minimal Residual method，QMR)，實現了大型稀疏線性方程組高效，精確的求 解。
[0018] 本發明與現有技術相比，有益效果在于:本發明針對實際大地電磁探測地下介質 中廣泛存在的各向異性問題，假設一個電性主軸垂直于層面方面，另外一個主軸平行于層 面方向，構造了各向異性介質模型的電導率張量模型。針對傳統的正演模擬方法依賴于網 格，物性參數復雜分布適應性差的缺點，提出用無網格徑向基點插值法模擬大地電磁的各 向異性問題，可以有效地克服單純多項式基點插值方法中存在的奇異性問題，形函數光滑 穩定，實現了電磁法高精度，自適應的數值模擬。
【附圖說明】
[0019] 圖1是本發明實施例提供的無網格法場節點、高斯點、求解域及其邊界，支持域與 背景網格示意圖；
[0020] 圖2是本發明實施例提供的各向異性介質的電性主軸與地下空間坐標系示意圖；
[0021] 圖3是本發明實施例提供的各向異性介質大地電磁無網格徑向基點插值法數值模 擬流程圖；
[0022]圖4是本發明實施例建立的層狀各向異性介質模型，其中旋轉角度為30° ;
[0023] 圖5是本發明實施例各項異性層狀模型擬解析解與無網格解的對比；
[0024] 圖6是本發明實施例建立的均勻介質中含有各向異性異常體模型；
[0025] 圖7是本發明實施例建立的均勻介質中含有各向異性異常體模型各個方向視電阻 率計算結果。
【具體實施方式】
[0026] 為了使本發明的目的、技術方案及優點更加清楚明白，以下結合實施例，對本發明 進行進一步詳細說明。應當理解，此處所描述的具體實施例僅僅用以解釋本發明，并不用于 限定本發明。
[0027] 參見圖1，本發明實施例提供的無網格法計算點1、求解域2,背景單元3,場節點4， 支持域5與求解域邊界6示意圖。
[0028] 求解域是整個邊值問題所要求解的區域，支持域是構造無網格法形函數所選擇的 區域，計算點是支持域的中心。在計算中，計算點常常選擇為高斯積分點。背景單元是便于 進行域積分而劃分的區域，背景單元并不依賴節點而存在。
[0029] 參見圖3,各向異性介質大地電磁無網格數值模擬方法流程圖，包括如下的步驟：
[0030] 1)從各向異性大地電磁邊值問題出發，構造等價泛函，推導對應無網格徑向基點 插值法的等價線性方程組；
[0031 ] 2)讀取當前模型參數，包括頻率參數，節點坐標，背景單元，支持域，形狀參數，極 化模式等；
[0032] 3)對當前極化模式進行判斷，若為TE極化模式，則計算區域包含空氣層，若為TM極 化模式，則不含空氣層；
[0033] 4)對所有背景網格進行循環，對背景網格的所有高斯積分點循環，搜索該高斯積 分點支持域內的有效節點，計算支持域內節點處的形函數，求取系數矩陣和右端項；
[0034] 5)加載本質邊界條件，求解線性方程組，得到各個節點的場值；
[0035] 6)由視電阻率計算公式，代入場值求取地面處各個方向的視電阻率和相位。
[0036] 當步驟4中計算系數矩陣和右端項結束后檢測高斯積分點循環是否結束，結束后 進一步檢測背景網格是否循環完畢，高斯積分點循環未結束時，返回步驟4對高斯積分點循 環，背景網格循環完畢后進入下一步，當背景網格未循環完畢時，則返回到步驟4對背景網 絡進行循環。
[0037]步驟1中各向異性大地電磁邊值問題 [0038]在介質中，頻率域電磁場方程是： VxE -iconII /
[0039] '' 、 ( 1 ) V x // = ^a-ico〇) E
[0040] 在地下介質空間中，對大地電磁探測所涉及的頻率而言，電導率〇>> e，因而也可 不考慮e的各向異性問題。電導率〇在各向同性介質中是標量，在各向異性介質中是張量。本 發明只考慮〇的各向異性問題。
[0041] 如圖2所示，建立如圖所示的坐標系，設平行層面的電導率為〇//，垂直層面的電導 率為〇丄，在層面坐標系x ' y ' z '中x '平行層面，y '垂直層面，z '平行走向，電導率張量為
(2)
[0043] 在地面坐標系中，取走向為z軸（與z'平行），x軸與z軸垂直，保持水平，y軸垂直向 上。在坐標系xyz中，電導率〇張量為
[0044] o=A〇，At (3)
[0045] 其中A為坐標變換張量
(4)
[0047]其中a是x'軸與地下空間坐標系x軸的夾角。
[0048]將（1)式按分量展開，并考慮到則將電導率張量代入，引入二維算子
經化簡可得各向異性大地電磁邊值問題如下： (5)
[0050]其中U指的是場函數Ez或者Hz，Q指的是求解域，AB是求解域的上邊界，AC和BD指的 是求解域左右邊界，CD指的是求解域下邊界，式中：
《代表圓頻率，y代表介 質的磁導率，在TE極化模式下，u = Ez，A = 〇//-i ? e :。在TM極化模式下，u = Hz，A=i ? y，T的表達式如下：
[0053 ]步驟4中支持域內徑向基點插值法建立形函數
[0054] 如圖1所示，在求解域Q內，徑向基點插值函數(RPIM)可表示為
[0056]式中Ri(X)為徑向基函數(RBF)，n為RBFs的個數，Pj(X)為空間坐標XT=(x，y)中的 單項式，m為多項式基函數的個數。&1和匕為待定常數。本發明選用的徑向基函數為復合2次 徑向基函數(MQ-RBF)，其表達式為：
[0057] Ri(x,y) = (ri2+(acdc)2)q(ac 0) (9)
[0058] 為確定式(8)中的ai和bj，需形成計算點X的支持域，其中包括n個場節點。使式(8) 滿足計算點X周圍的n個節點值以確定系數 &1和匕，這將