基于非結(jié)構(gòu)化網(wǎng)格的大規(guī)模磁測(cè)數(shù)據(jù)正則化反演研究

劉陽(yáng)，張志勇

（東華理工大學(xué) 地球物理與測(cè)控技術(shù)學(xué)院，江西 南昌 330013）

磁法勘探已廣泛應(yīng)用于鈾礦勘探、深部地質(zhì)構(gòu)造圈定、地?zé)嵴{(diào)查等領(lǐng)域［1-4］，大規(guī)模鈾礦地質(zhì)勘探亟需準(zhǔn)確、高效的磁法解譯。位場(chǎng)正演是磁法勘探解釋任務(wù)的基礎(chǔ)，其計(jì)算成本依舊是阻礙三維反演發(fā)展的一個(gè)嚴(yán)重問(wèn)題［5］。位場(chǎng)正演模擬主要有空間域和頻率域兩類(lèi)方法，其中空間域方法在處理起伏地表與非規(guī)則離散模型等方面具有較大優(yōu)勢(shì)。空間域方法又可分為解析法和數(shù)值法。通過(guò)閉合積分公式建模的空間域方法已有眾多研究成果［6-7］，簡(jiǎn)單矩形棱柱積分公式難以擬合復(fù)雜地質(zhì)體，均勻多面體剖分算法逐漸成為位場(chǎng)正演建模的主流［8］。然而，三維復(fù)雜地質(zhì)模型計(jì)算公式復(fù)雜，適應(yīng)性較差、計(jì)算量較大［9］，以解析法正演開(kāi)展大規(guī)模離散反演時(shí)需要大量積分運(yùn)算和內(nèi)存消耗，計(jì)算成本巨大，因此，位場(chǎng)解析法建模的應(yīng)用受到了限制。

位場(chǎng)數(shù)值法建模主要有有限差分法［10］、有限體積法［11-12］、積分方程法［13］、有限單元法等［14］，與解析法相比，數(shù)值法計(jì)算精度略低，但計(jì)算時(shí)間和內(nèi)存消耗方面具有一定優(yōu)勢(shì)。其中有限差分法通常采用結(jié)構(gòu)化網(wǎng)格，但在物性參數(shù)復(fù)雜分布或場(chǎng)域分布幾何特征不規(guī)則時(shí)適應(yīng)性差，存在一定局限性。有限單元法既可采用結(jié)構(gòu)化網(wǎng)格，也可采用非結(jié)構(gòu)化網(wǎng)格。與結(jié)構(gòu)化網(wǎng)格相比，非結(jié)構(gòu)化網(wǎng)格在要求較高精度或可接受較低精度的區(qū)域，網(wǎng)格的局部細(xì)化或粗化可以在不影響這些區(qū)域以外的網(wǎng)格的情況下執(zhí)行［12］。因此，使用非結(jié)構(gòu)化網(wǎng)格的數(shù)值格式會(huì)有更高的精度、更低的計(jì)算時(shí)間和內(nèi)存占用。有限單元法以其適用性強(qiáng)、可以模擬任意復(fù)雜地形及地質(zhì)體而受到廣泛運(yùn)用［15-17］，如地震［18］、直流電阻率［15］、自然電位［19］、大地電磁［20］正反演中均得到應(yīng)用。前人基于有限單元法開(kāi)展了大量位場(chǎng)正反演的研究，Zhang等［21］將遺傳算法與有限單元法相結(jié)合實(shí)現(xiàn)了一種替代方法反演重力異常數(shù)據(jù)，重建地下三維密度結(jié)構(gòu)；Cai 等［14］通過(guò)一種快速的有限單元法來(lái)解決復(fù)雜地質(zhì)條件下引力位的邊值問(wèn)題；Maag 等［22］由重力位的偏微分方程出發(fā)推導(dǎo)出基于有限單元法的正演和反演問(wèn)題；Jahandari 和Farquharson［12］采用有限單元法和有限體積法與解析法對(duì)比實(shí)現(xiàn)了重力位正反演；但當(dāng)前文獻(xiàn)少有以非結(jié)構(gòu)化有限單元法開(kāi)展磁場(chǎng)正反演的研究。

大規(guī)模磁測(cè)數(shù)據(jù)反演中，靈敏度矩陣的儲(chǔ)存與計(jì)算一定程度上制約著計(jì)算規(guī)模的大小。通過(guò)將靈敏度矩陣及其轉(zhuǎn)置與向量乘積轉(zhuǎn)換為正演過(guò)程，進(jìn)行有限內(nèi)存最優(yōu)化求解的反演策略已應(yīng)用到地球物理反演問(wèn)題中，其中算法包括共軛梯度法［23］、有限內(nèi)存擬牛頓法［24-25］、高斯-牛頓法［26］等，與其他最優(yōu)化算法相比，高斯-牛頓法在保證快速收斂的同時(shí)，有效減少了計(jì)算量，因此，高斯-牛頓法適用于大規(guī)模磁測(cè)數(shù)據(jù)反演。

針對(duì)以解析法正演為依托的磁異常反演及規(guī)則網(wǎng)格數(shù)值解對(duì)復(fù)雜地質(zhì)體擬合程度不足、計(jì)算成本大等缺點(diǎn)，本文采用非結(jié)構(gòu)化網(wǎng)格有限單元法進(jìn)行磁異常正演，通過(guò)靈敏度矩陣隱式存儲(chǔ)的高斯-牛頓方法求解反演目標(biāo)函數(shù)實(shí)現(xiàn)磁測(cè)數(shù)據(jù)正則化反演。研究工作首先采用局部網(wǎng)格加密技術(shù)和最小二乘空間梯度計(jì)算方法提高有限單元正演計(jì)算精度，通過(guò)與解析法對(duì)比，對(duì)非結(jié)構(gòu)化網(wǎng)格有限單元法正演精度的影響因素進(jìn)行了試算討論；在此基礎(chǔ)之上，開(kāi)展了理論模型與鈾礦區(qū)域?qū)崪y(cè)數(shù)據(jù)的反演，驗(yàn)證了本文提出的三維磁測(cè)數(shù)據(jù)反演算法的可行性。

1 三維有限單元正演

1.1 正演理論

在假設(shè)無(wú)傳導(dǎo)電流的空間域，磁位基本方程可表示為：

式中：μ0為真空磁導(dǎo)，H/m；U為磁位；κ為磁化率；?為哈密爾頓算子；Mr為剩余磁化強(qiáng)度矢量，A/m。當(dāng)在兩種介質(zhì)存在的空間域中，采用有限單元法求解方程式（1），考慮到兩種介質(zhì)分界面磁位連續(xù)且無(wú)散場(chǎng)的特點(diǎn)，得到公式：

式中：n為分界面上的單位正法線(xiàn)方向；μ1和μ2為邊界兩種介質(zhì)中的磁導(dǎo)率；U1和U2為邊界兩種介質(zhì)中磁位；Mr1和Mr2為邊界兩側(cè)介質(zhì)中剩余磁化強(qiáng)度矢量，Mr1，n和Mr2，n是Mr1和Mr2在 外法向的投影，則有第一類(lèi)邊界條件：

式中：T為總磁場(chǎng)強(qiáng)度矢量，T；r為坐標(biāo)原點(diǎn)到無(wú)窮遠(yuǎn)邊界距離矢量。第二類(lèi)邊界條件可表示為：

磁場(chǎng)邊值問(wèn)題與求解下式變分問(wèn)題等價(jià)：

式中：Tn是磁場(chǎng)T在邊界外法線(xiàn)方向上的投影，Ω 為研究區(qū)域。對(duì)計(jì)算區(qū)域進(jìn)行非結(jié)構(gòu)化四面體離散，在每個(gè)剖分單元采用線(xiàn)性函數(shù)插值，對(duì)各單元泛函求和，得到下列線(xiàn)性系統(tǒng)：

式中：K為剛度矩陣（對(duì)稱(chēng)稀疏矩陣）；U為節(jié)點(diǎn)磁位；P為邊界條件形成源項(xiàng)。求解聯(lián)立方程組可得到各節(jié)點(diǎn)磁位U，通過(guò)數(shù)值微商可得到各磁場(chǎng)分量，但由于非結(jié)構(gòu)化網(wǎng)格中網(wǎng)格大小不一致，基于線(xiàn)性插值計(jì)算磁位，對(duì)磁位進(jìn)行求導(dǎo)磁場(chǎng)為常數(shù)，使得在計(jì)算磁場(chǎng)時(shí)的精度和穩(wěn)定性存在問(wèn)題。本文擬采用與觀測(cè)點(diǎn)共點(diǎn)的所有四面體進(jìn)行空間梯度的最小二乘計(jì)算。假設(shè)定義觀測(cè)點(diǎn)坐標(biāo)為(x0，y0，z0)，相應(yīng)磁位為U，與觀測(cè)點(diǎn)共點(diǎn)的所有四面體n 個(gè)頂點(diǎn)坐標(biāo)定義為(xi，yi，zi)，相應(yīng)磁位為Ui，其中i=1，2…n，則有：

節(jié)點(diǎn)組成方程組為：

可將式（9）表示為：

令E=(U-AL)T(U-AL)，對(duì)L 進(jìn)行求導(dǎo)可得：

令上式為0，得L=(AT A)-1ATU，求解方程可得到坐標(biāo)點(diǎn)處梯度為：

進(jìn)而可得到觀測(cè)點(diǎn)處磁場(chǎng)強(qiáng)度。

1.2 正演精度影響因素分析

正演精度直接影響反演的可靠性，對(duì)于非結(jié)構(gòu)化有限單元正演，其精度受多種因素影響。本文從觀測(cè)網(wǎng)格密度、邊界節(jié)點(diǎn)數(shù)兩方面進(jìn)行討論，并與均勻多面體的閉合積分公式［27］計(jì)算出的總磁異常進(jìn)行精度對(duì)比。設(shè)計(jì)模型如圖1所示，假設(shè)在地下空間磁化率為0.006 SI 的背景中，設(shè)置一個(gè)磁化率為0.1 SI 的異常體，異常體頂板埋深為500 m，結(jié)構(gòu)為1 000 m×1 000 m×1 000 m，地磁場(chǎng)強(qiáng)度在X、Y、Z 三個(gè)方向的分量分別為（30 000 nT，0 nT，40 000 nT）。選取觀測(cè)區(qū)任意部分區(qū)域測(cè)點(diǎn)驗(yàn)證有限單元法的計(jì)算精度，計(jì)算區(qū)域大小為250 m×250 m。

圖1 磁化率模型Fig.1 Magnetic susceptibility model

1.2.1 觀測(cè)網(wǎng)格影響

在有限單元計(jì)算中，以觀測(cè)點(diǎn)作為剖分網(wǎng)格中心點(diǎn)進(jìn)行剖分，因此，有限單元計(jì)算精度依賴(lài)于剖分網(wǎng)格到中心點(diǎn)的距離。對(duì)不同觀測(cè)網(wǎng)格產(chǎn)生總磁異常進(jìn)行計(jì)算，邊界節(jié)點(diǎn)數(shù)為38，計(jì)算結(jié)果如圖2 所示。不同觀測(cè)網(wǎng)格的剖分網(wǎng)格屬性及各類(lèi)誤差信息如表1 所示。隨著觀測(cè)網(wǎng)格變密，平均絕對(duì)誤差從3.67 nT降至1.60 nT，最大絕對(duì)誤差從14.08 nT 下降至5.56 nT，平均相對(duì)誤差從0.047 下降至0.018，計(jì)算總磁異常精度明顯提高。

表1 地表不同觀測(cè)密度所產(chǎn)生的網(wǎng)格參數(shù)及各類(lèi)誤差信息表Table 1 Grid parameters and error information of different observing density

圖2 不同觀測(cè)網(wǎng)格下的總磁異常等值線(xiàn)及絕對(duì)誤差圖Fig.2 Contour of total magnetic anomaly and absolute error under different observation grids

1.2.2 邊界節(jié)點(diǎn)數(shù)影響

為分析邊界節(jié)點(diǎn)數(shù)對(duì)有限單元正演模擬計(jì)算的影響，選取邊界節(jié)點(diǎn)［15，25，35，45］進(jìn)行試算，觀測(cè)網(wǎng)格為50 m×50 m，以1.2 倍指數(shù)擴(kuò)邊得到圖3 和表2。計(jì)算結(jié)果顯示，隨著邊界節(jié)點(diǎn)數(shù)的增加，平均絕對(duì)誤差從80.15 nT 下降至3.77 nT，最大絕對(duì)誤差從95.33 nT 下降至12.31 nT，平均相對(duì)誤差從0.460 下降至0.046，最大相對(duì)誤差減小，計(jì)算精度明顯提高，而過(guò)少邊界節(jié)點(diǎn)計(jì)算總磁異常存在較大的誤差。

表2 不同邊界節(jié)點(diǎn)數(shù)計(jì)算產(chǎn)生網(wǎng)格參數(shù)及各類(lèi)誤差信息表Table 2 Grid parameters and error information calculated with different boundary nodes

圖3 不同邊界節(jié)點(diǎn)數(shù)計(jì)算總磁異常等值線(xiàn)及絕對(duì)誤差圖Fig.3 Contour of total magnetic anomaly and absolute error calculated with different boundary nodes

1.2.3 局部網(wǎng)格密度影響

觀測(cè)網(wǎng)格密度和邊界節(jié)點(diǎn)數(shù)對(duì)于有限單元計(jì)算精度均有影響，而觀測(cè)網(wǎng)格變密導(dǎo)致計(jì)算效率的下降，此外，觀測(cè)網(wǎng)格密度達(dá)到一定程度，繼續(xù)加密對(duì)于計(jì)算精度有較小的影響［15］。因此，本文采用網(wǎng)格加密的方式提高磁異常計(jì)算精度，網(wǎng)格加密的方法主要有局部體積加密方法（Local Volume-Refinement，簡(jiǎn)稱(chēng)LVR）、全局體積加密方法（Global Volume-Refinement，簡(jiǎn)稱(chēng)GVR）和局部測(cè)點(diǎn)加密方法（Local Node-Refinement，簡(jiǎn)稱(chēng)LNR）。其中GVR 法加密基于結(jié)構(gòu)化網(wǎng)格會(huì)導(dǎo)致過(guò)多的冗余加密節(jié)點(diǎn)，精度提高不明顯；LVR 法是在加密需求區(qū)域無(wú)應(yīng)有的加密，反而在其他區(qū)域加入過(guò)多節(jié)點(diǎn)，導(dǎo)致計(jì)算量增大；LNR法相對(duì)于LVR 和GVR 法可以減少內(nèi)存消耗和計(jì)算時(shí)間，因此，本文采用LNR 網(wǎng)格加密方式。

局部測(cè)點(diǎn)網(wǎng)格細(xì)化策略是基于觀測(cè)點(diǎn)位置下方插入節(jié)點(diǎn)，觀測(cè)網(wǎng)格為50 m×50 m，擴(kuò)邊節(jié)點(diǎn)數(shù)為45，得到計(jì)算結(jié)果如圖4 和表3 所示。局部網(wǎng)格細(xì)化可提高正演計(jì)算精度，插入節(jié)點(diǎn)深度較淺時(shí)能獲得更高精度的異常分辨率，當(dāng)插入節(jié)點(diǎn)后剖分網(wǎng)格數(shù)量增加明顯，而插入節(jié)點(diǎn)接近觀測(cè)點(diǎn)時(shí)剖分網(wǎng)格數(shù)量明顯多于插入節(jié)點(diǎn)較深的情況，當(dāng)在觀測(cè)點(diǎn)下60 m、20 m、1 m 插入節(jié)點(diǎn)時(shí)平均絕對(duì)誤差從2.47 nT減小至1.11 nT，平均相對(duì)誤差從0.028 下降至0.011，相比同一觀測(cè)密度計(jì)算，局部網(wǎng)格加密后磁異常精度和穩(wěn)定性明顯提高。

表3 不同深度插入節(jié)點(diǎn)產(chǎn)生網(wǎng)格參數(shù)及各類(lèi)誤差信息表Table 3 Mesh parameters and various error information calculated by inserting nodes at different depths

圖4 不同深度插入節(jié)點(diǎn)計(jì)算磁異常等值線(xiàn)及絕對(duì)誤差圖Fig.4 Contour of magnetic anomaly and absolute error calculated by inserting nodes at different depths

以上試算結(jié)果表明，通過(guò)提高觀測(cè)網(wǎng)格密度、增多邊界節(jié)點(diǎn)數(shù)及采用局部網(wǎng)格細(xì)化技術(shù)等可實(shí)現(xiàn)非結(jié)構(gòu)化網(wǎng)格有限單元法計(jì)算磁異常精確解。

1.3 有限單元計(jì)算效率分析

為說(shuō)明非結(jié)構(gòu)化網(wǎng)格有限單元法計(jì)算的低耗存和高效率的特點(diǎn)，本文測(cè)試平臺(tái)為Intel（R）Core（TM）i7-9700F CPU@3.00 GHz，RAM 32.0 GB。在同樣的觀測(cè)網(wǎng)格及相同剖分網(wǎng)格參數(shù)下，進(jìn)行有限單元法與解析法計(jì)算內(nèi)存和時(shí)長(zhǎng)的統(tǒng)計(jì)（表4 和圖5）。如表4 和圖5 所示，與解析算法相比，同樣的計(jì)算條件和剖分網(wǎng)格參數(shù)下，有限單元法計(jì)算時(shí)長(zhǎng)與內(nèi)存需求明顯優(yōu)于解析算法。

表4 有限單元法與解析法計(jì)算效率對(duì)比Table 4 Comparison of computational efficiency between finite element method and analytical method

圖5 計(jì)算內(nèi)存與時(shí)長(zhǎng)曲線(xiàn)對(duì)比圖Fig.5 Contrast of calculation memory versus duration curve between the method of finite element and analysis

2 大規(guī)模磁測(cè)數(shù)據(jù)反演算法

2.1 正則化目標(biāo)函數(shù)

正則化反演目標(biāo)函數(shù)一般形式可表示為［28］：

式中：?d(m)為數(shù)據(jù)誤差函數(shù)；?m(m)為模型穩(wěn)定函數(shù)，本次研究工作采用最小結(jié)構(gòu)穩(wěn)定函數(shù)［29］；α為正則化因子，是平衡數(shù)據(jù)誤差函數(shù)和模型穩(wěn)定函數(shù)之間的權(quán)重系數(shù)，采用L-Curve 方法選?。?0］；dobs為觀測(cè)數(shù)據(jù)向量；m為離散模型向量；mapr為先驗(yàn)?zāi)Ｐ蛥?shù)向量；A為正演算子；Wd為數(shù)據(jù)協(xié)方差矩陣；Wm為模型權(quán)重矩陣。

2.2 目標(biāo)函數(shù)優(yōu)化

本次研究工作采用高斯-牛頓方法求解目標(biāo)函數(shù)：

式（13）數(shù)據(jù)誤差項(xiàng)和模型穩(wěn)定項(xiàng)展開(kāi)可得：

將上述方程對(duì)Δm(n)求導(dǎo)，可得高斯-牛頓方程：

采用J表示靈敏度矩陣，J中任意元素表示為：

式中：Hi表示第i 個(gè)觀測(cè)數(shù)據(jù)；m表示反演模型參數(shù)；mj為第j 個(gè)反演模型參數(shù)。將其與有限單元線(xiàn)性方程式（7）結(jié)合，則：

將上式代入式（19），可得：

將靈敏度矩陣及其轉(zhuǎn)置與任意矩陣的乘積轉(zhuǎn)換成線(xiàn)性方程的求解，并且該方程與有限單元具有相同的矩陣：

由式（25）和（26）所示，求取靈敏度矩陣（或轉(zhuǎn)置）和任意矩陣的乘積可包含在靈敏度求解過(guò)程中，因此，實(shí)現(xiàn)了靈敏度矩陣的隱式儲(chǔ)存，減少了內(nèi)存消耗，提高了計(jì)算效率。

3 模型試算分析

3.1 單一模型

首先采用單一模型驗(yàn)證本文提出算法的可行性，單一模型和復(fù)雜模型正演數(shù)據(jù)中均加入5%隨機(jī)噪聲。模型設(shè)置如圖1，采用總磁異常數(shù)據(jù)進(jìn)行反演，觀測(cè)網(wǎng)格為50 m×50 m，測(cè)點(diǎn)數(shù)為4 900 個(gè)，初始正則化因子為1.0×10-5，迭代75 次。如圖6 所示，采用本文算法對(duì)單一模型反演取得了良好的效果，其中白色矩形框?yàn)槟Ｐ臀恢茫囱萁Y(jié)果磁化率分布與真實(shí)模型位置有較好對(duì)應(yīng)。從RMS 誤差曲線(xiàn)可知（圖7），隨著迭代的進(jìn)行，RMS 總體趨勢(shì)穩(wěn)定下降趨近于1.0，說(shuō)明反演穩(wěn)定收斂；反演中通過(guò)將靈敏度矩陣及其轉(zhuǎn)置與任意向量的乘積轉(zhuǎn)換成正演計(jì)算，內(nèi)存需求減少，計(jì)算效率明顯提高。

圖6 單一模型反演結(jié)果Fig.6 Inversion results of single model

圖7 單一模型反演RMS 誤差曲線(xiàn)Fig.7 RMS error curve of inversion result by single model

3.2 復(fù)雜模型

為進(jìn)一步測(cè)試本文提出的反演方法對(duì)復(fù)雜模型的反演效果，設(shè)置復(fù)雜磁化率模型如圖8 所示。假設(shè)在地下空間磁化率為0.006 SI 的背景中，設(shè)置4 個(gè)磁化率均為0.1 SI 的矩形異常體，異常體結(jié)構(gòu)及屬性見(jiàn)表5，觀測(cè)網(wǎng)格為50 m×50 m，測(cè)點(diǎn)數(shù)共計(jì)4 900 個(gè)，地磁場(chǎng)強(qiáng)度與單一模型等同。對(duì)于復(fù)雜模型進(jìn)行磁化率反演取得同樣效果（圖9），反演磁化率分布位置與設(shè)定的4 個(gè)異常體模型位置均有對(duì)應(yīng)，且反演磁化率橫向分辨率有較好體現(xiàn)（圖9a）。從復(fù)雜模型反演RMS 誤差曲線(xiàn)可以發(fā)現(xiàn)（圖10），RMS 隨迭代次數(shù)逐漸下降并最終趨近于1.0，說(shuō)明反演穩(wěn)定收斂；相比解析法，以有限單元法開(kāi)展復(fù)雜模型離散反演計(jì)算效率明顯提高，試算結(jié)果表明，本文提出反演算法對(duì)于復(fù)雜模型同樣具有可行性。

表5 異常體屬性Table 5 Attributes of anomaly body

圖8 設(shè)置磁化率模型Fig.8 Setting of the magnetic susceptibility model

圖9 復(fù)雜模型反演結(jié)果Fig.9 Inversion results of complex model

圖10 復(fù)雜模型反演RMS 誤差曲線(xiàn)Fig.10 RMS error curve of inversion result by complex model

3.3 實(shí)測(cè)數(shù)據(jù)反演

磁法在鈾礦地質(zhì)勘查中應(yīng)用廣泛，為驗(yàn)證本文提出算法處理鈾礦區(qū)域采集實(shí)際數(shù)據(jù)的有效性，選取某鈾礦研究區(qū)帶地形磁測(cè)數(shù)據(jù)進(jìn)行反演。研究區(qū)巖性以花崗巖為主，勘探目標(biāo)為輝綠巖，觀測(cè)網(wǎng)格為25 m×25 m，測(cè)點(diǎn)數(shù)共計(jì)3 025 個(gè)，研究區(qū)反演最高處為15 m，磁場(chǎng)強(qiáng)度為當(dāng)?shù)氐卮艌?chǎng)，磁異常為輝綠巖引起。初始正則化因子為1.0×10-8，共計(jì)迭代28 次。圖11 a 為實(shí)測(cè)數(shù)據(jù)反演網(wǎng)格結(jié)果，反演節(jié)點(diǎn)數(shù)為62 841 個(gè)，四面體單元數(shù)為390 828 個(gè)，三角單元數(shù)為783 079 個(gè)，邊數(shù)為455 091 個(gè)；圖11b為反演結(jié)果橫切圖，切片深度分別為［10 m，-40 m，-70 m，-110 m］。

圖11 實(shí)測(cè)數(shù)據(jù)反演結(jié)果Fig.11 Inversion results of measured data

如圖11 所示，應(yīng)用本文提出算法對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行處理得到較好的磁化率分布特征，模型邊界分辨率高，異常特征明顯，基于非結(jié)構(gòu)網(wǎng)格其適用性強(qiáng)，具有模擬復(fù)雜結(jié)構(gòu)的能力，對(duì)實(shí)測(cè)數(shù)據(jù)地形有較好的擬合。從實(shí)測(cè)數(shù)據(jù)反演RMS 誤差曲線(xiàn)可以看出（圖12），隨著迭代次數(shù)RMS 整體呈下降趨勢(shì)，說(shuō)明實(shí)測(cè)數(shù)據(jù)反演同樣穩(wěn)定收斂?；谝陨显囁阏f(shuō)明應(yīng)用本文提出算法進(jìn)行大規(guī)模帶地形磁測(cè)數(shù)據(jù)反演具有有效性，為大區(qū)域鈾礦地質(zhì)勘查提供一定技術(shù)支持。

圖12 實(shí)測(cè)數(shù)據(jù)反演RMS 誤差曲線(xiàn)Fig.12 RMS error curve of inversion result from measured data

4 結(jié)論

本文開(kāi)展了基于非結(jié)構(gòu)化網(wǎng)格的大規(guī)模磁測(cè)數(shù)據(jù)三維正反演研究，分析了觀測(cè)網(wǎng)格密度、邊界節(jié)點(diǎn)數(shù)及局部網(wǎng)格細(xì)化等因素對(duì)正演精度的影響，探討了利用靈敏度矩陣隱式存儲(chǔ)的方法實(shí)現(xiàn)大規(guī)模磁測(cè)數(shù)據(jù)反演的算法，得到以下結(jié)論：

1）磁異常的非結(jié)構(gòu)化網(wǎng)格的有限單元正演在內(nèi)存消耗和計(jì)算效率兩方面均優(yōu)于解析算法。

2）非結(jié)構(gòu)化網(wǎng)格有限單元正演的網(wǎng)格離散參數(shù)對(duì)計(jì)算精度有直接影響。提高觀測(cè)網(wǎng)格密度、增多邊界節(jié)點(diǎn)數(shù)并采用局部網(wǎng)格細(xì)化技術(shù)，可提高磁測(cè)數(shù)據(jù)正演的精度和穩(wěn)定性。

3）非結(jié)構(gòu)網(wǎng)格有限單元正演與靈敏度矩陣隱式存儲(chǔ)的高斯-牛頓法反演相結(jié)合，有效減小了反演的內(nèi)存消耗，降低了計(jì)算成本，適用于大規(guī)模三維磁測(cè)數(shù)據(jù)反演。