二維直流電阻率與音頻大地電磁自適應(yīng)漸進(jìn)聯(lián)合反演

李 曼，于 鵬，張志勇，張寶松

（1.同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院，上海200092；2.東華理工大學(xué)地球物理與測控技術(shù)學(xué)院，江西南昌330013；3.中國地質(zhì)調(diào)查局南京地質(zhì)調(diào)查中心，江蘇南京210016）

聯(lián)合反演可有效降低反演的多解性，提高反演結(jié)果的精度和可靠性［1-2］。直流電阻率（DC）與大地電磁（MT）聯(lián)合反演研究始于20世紀(jì)70年代，一維研究表明，聯(lián)合反演可避免單一方法反演的不確定性［3］；Constable等［4］采用平滑約束的直流測深與MT一維歐克姆（OCCAM）反演研究表明，MT、DC數(shù)據(jù)誤差不一致會(huì)造成兩個(gè)數(shù)據(jù)集擬合程度的不同，同時(shí)實(shí)際模型的非一維特征對反演存在影響。在二維方面，基于矩形剖分網(wǎng)格、平滑模型約束、互換定理計(jì)算偏導(dǎo)數(shù)矩陣的MT與DC聯(lián)合反演實(shí)踐表明，聯(lián)合反演可以減少單獨(dú)反演的不確定性［5-6］。類似MT方法，通過觀測正交的電場與磁場分量，計(jì)算阻抗并定義視電阻率的射頻電磁法即RMT方法（radiomagnetotelluric，頻段10 kHz～1 MHz），由于其勘探深度與直流電阻率有很大的重疊，推動(dòng)了聯(lián)合反演工作。研究表明，聯(lián)合反演RMT與DC數(shù)據(jù)可解決靜態(tài)效應(yīng)對RMT的影響［7］。單獨(dú)反演與聯(lián)合反演的對比研究表明，DC對高阻與低阻地質(zhì)體均有較好響應(yīng)，而RMT對低阻體更靈敏，通過合理選擇兩個(gè)數(shù)據(jù)集在反演中的權(quán)重可以得到更佳的反演結(jié)果［8］。實(shí)踐工作表明，二維聯(lián)合反演可以得到比三維張量反演更精細(xì)的地下結(jié)構(gòu)［9］。采用相同模型約束條件下的概率反演與確定性反演的效果相當(dāng)，聯(lián)合反演改善了RMT反演效果［10］。

鑒于DC與MT、DC與RMT聯(lián)合反演可有效改善反演效果，本文進(jìn)一步開展非結(jié)構(gòu)網(wǎng)格剖分條件下DC與音頻大地電磁（AMT）二維聯(lián)合反演的研究工作。開發(fā)了以模型靈敏度信息為依據(jù)的反演網(wǎng)格優(yōu)化技術(shù)，以生成高質(zhì)量反演網(wǎng)格；構(gòu)建由粗網(wǎng)格到細(xì)網(wǎng)格的漸進(jìn)反演策略，減少反演對穩(wěn)定因子的依賴，從而簡化正則化因子計(jì)算搜索的算法；通過最小二乘方法求解模型空間變化趨勢方程系數(shù)，進(jìn)而求解非結(jié)構(gòu)三角網(wǎng)格的二階梯度算子，構(gòu)建非結(jié)構(gòu)化網(wǎng)格條件下最小結(jié)構(gòu)穩(wěn)定因子；采用高斯-牛頓法優(yōu)化求解正則化目標(biāo)函數(shù)，利用雙共軛梯度穩(wěn)定算法（BICGSTAB）求解高斯-牛頓方程確保反演穩(wěn)定，減少內(nèi)存需求；最后，通過合成數(shù)據(jù)與實(shí)測數(shù)據(jù)的反演計(jì)算，驗(yàn)證了上述方法的有效性。

1 DC與AMT正則化聯(lián)合反演

為更好地?cái)M合地表與地下地質(zhì)單元，實(shí)現(xiàn)野外工作點(diǎn)距差別大、勘探深度與分辨率不同的AMT與DC兩種方法的高效、高精度正演，本文采用非結(jié)構(gòu)化三角網(wǎng)格進(jìn)行區(qū)域剖分；應(yīng)用有限單元法進(jìn)行DC與AMT正演［11-14］；采用基于最小結(jié)構(gòu)模型的正則化方法構(gòu)造非結(jié)構(gòu)化三角網(wǎng)格條件下的模型二階平滑計(jì)算方法，采用高斯-牛頓最優(yōu)化方法求解反演目標(biāo)函數(shù)，并通過雙共軛梯度穩(wěn)定算法（BICGSTAB）不完全求解高斯-牛頓方程，以確保反演過程穩(wěn)定，同時(shí)減少反演對內(nèi)存的需求。

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

地球物理反問題通常為病態(tài)問題，正則化是穩(wěn)定求解反問題的有效手段［15-16］。正則化反演的目標(biāo)函數(shù)為

式中：m為反演參數(shù)向量；dobs為反演數(shù)據(jù)向量；μ為正則化因子；φd(m，dobs)為數(shù)據(jù)誤差函數(shù)；φm(m)為模型誤差函數(shù)。

DC與AMT聯(lián)合反演的變量為介質(zhì)電阻率，采用以下變換函數(shù)［17-18］以確保反演電阻率在實(shí)際巖、礦石的取值范圍內(nèi)：

式中：mi為反演空間參數(shù)；ρi為反演單元的電阻率值；ρmin、ρmax分別為研究區(qū)電阻率的下限和上限。

反演數(shù)據(jù)由AMT與DC數(shù)據(jù)組成，dobs={ρTM，a，ρTE，a，φTM，a，φTE，a，ρDC，s}，其中 ，ρTM，a、φTM，a分別為音頻大地電磁TM模式視電阻率與相位；ρTE，a、φTE，a分別為TE模式視電阻率與相位；ρDC，s為直流電阻率法觀測的視電阻率。數(shù)據(jù)誤差函數(shù)可表示為如下統(tǒng)一形式［16］：

式中：G(m)為正演函數(shù)，分別由AMT與DC兩種方法觀測數(shù)據(jù)組成；Wd為數(shù)據(jù)方差矩陣，

其中，ε為確保分母不為零的一個(gè)正實(shí)數(shù)；χi為數(shù)據(jù)方差；αi為數(shù)據(jù)權(quán)重系數(shù)，其作用是平衡AMT與DC數(shù)據(jù)在反演中的權(quán)重，以避免由于兩種數(shù)據(jù)間誤差范圍的不同造成某一數(shù)據(jù)集的權(quán)重過大，而另一數(shù)據(jù)不起作用。在數(shù)據(jù)誤差為高斯噪聲的假設(shè)下，反演后數(shù)據(jù)誤差函數(shù)期望值為NMT+NDC，NMT為AMT參與反演的數(shù)據(jù)個(gè)數(shù)、NDC為DC參與反演的數(shù)據(jù)個(gè)數(shù)。如果希望AMT與DC兩個(gè)數(shù)據(jù)集在反演中取相同比重，數(shù)據(jù)權(quán)重系數(shù)可表示為

先驗(yàn)信息通常根據(jù)模型空間分布特征進(jìn)行構(gòu)造，根據(jù)其反演后模型的特點(diǎn)可大致分為平滑與陡變兩種約束。其中平滑模型一般利用模型參數(shù)空間梯度進(jìn)行計(jì)算。支持陡變模型的約束方法有多種，如最大變化量支持模型、最小支持模型、最小梯度支持模型等；另外，通過模型誤差的L1與L2范數(shù)形式可以實(shí)現(xiàn)模型分塊與平滑形式的反演約束。為確保反演的穩(wěn)定性，本文采用最小結(jié)構(gòu)模型約束。

式中：βs，βx，βy為模型誤差三部分之間的比例系數(shù)；mapr為先驗(yàn)?zāi)Ｐ?；模型誤差可統(tǒng)一表示為

式中：Wm為統(tǒng)一定義形式的模型誤差矩陣。

式（1）中μ為平衡數(shù)據(jù)誤差與模型誤差的系數(shù)（正則化因子），可采用廣義交叉檢驗(yàn)（GCV）、“L－Curve”方法、貝葉斯方法、經(jīng)驗(yàn)選擇方法等進(jìn)行計(jì)算。正則化因子的選擇關(guān)系反演的成敗，若采用GCV、LCurve、貝葉斯方法等選取正則化因子，計(jì)算量較大，嚴(yán)重影響反演計(jì)算效率。本文研究了漸進(jìn)反演技術(shù)，在反演過程中盡可能減少對穩(wěn)定因子的依賴，因此采用簡單的經(jīng)驗(yàn)選擇方法即可實(shí)現(xiàn)μ的選擇。研究工作所采用的經(jīng)驗(yàn)選擇方法以數(shù)據(jù)誤差、數(shù)據(jù)誤差期望值等信息為依據(jù)，進(jìn)行正則化因子的動(dòng)態(tài)調(diào)整，基本策略如下：首次反演迭代設(shè)置正則化因子初值，設(shè)定數(shù)據(jù)誤差的期望值；當(dāng)?shù)螖?shù)大于1時(shí)，如數(shù)據(jù)誤差大于其期望值，按下式進(jìn)行調(diào)整：

式中：γ＞1.0為調(diào)整系數(shù)；Du、Dd為數(shù)據(jù)誤差下降速度的限制參數(shù)，當(dāng)數(shù)據(jù)誤差下降過快，增加正則化因子，若下降過慢則減少正則化因子，本次研究工作取Du=2.0，Dd=1.3。如數(shù)據(jù)誤差小于期望值，則按比例增大正則化因子；另外，考慮到正演計(jì)算誤差，在漸進(jìn)反演初期采用較大的數(shù)據(jù)誤差期望，隨反演網(wǎng)格的細(xì)化減小數(shù)據(jù)誤差期望。

1.2 非結(jié)構(gòu)化三角網(wǎng)格模型梯度計(jì)算

矩形剖分時(shí)，模型粗糙度一般通過水平與垂直方向差分計(jì)算；而非結(jié)構(gòu)三角網(wǎng)格，由于單元邊的法向方向一般不在水平與垂直方向且各單元之間也不一致，所以計(jì)算更為復(fù)雜。根據(jù)單元的空間梯度與方向梯度存在線性關(guān)系，Lelièvre等［19］給出了非結(jié)構(gòu)網(wǎng)格一階梯度的計(jì)算方法；Key［18］采用了與計(jì)算單元存在共用頂點(diǎn)的三角單元進(jìn)行粗糙度計(jì)算的方法，并在單元中心距離的計(jì)算中引入方向懲罰因子，實(shí)現(xiàn)水平與垂直方向不同強(qiáng)度的約束；?zyildirim等［20］采用了與計(jì)算單元具有共用邊的三角單元，并通過邊長與周長比值進(jìn)行差分定義的粗糙度計(jì)算方法。由于Lelièvre算法意義明確，可將其擴(kuò)展以實(shí)現(xiàn)二階梯度計(jì)算。取第i個(gè)三角單元中心坐標(biāo)為(xi，yi)，介質(zhì)物性值為ωi，如圖1所示。

圖1 非結(jié)構(gòu)三角網(wǎng)格梯度Fig.1 Gradient calculation of unstructured triangular mesh

假設(shè)物性分布為空間坐標(biāo)的二次函數(shù)，有

通過與計(jì)算單元具有共用頂點(diǎn)的三角形組成的集合進(jìn)行二階梯度計(jì)算，用最小二乘求解式（9）的系數(shù)為

式中：q=(a，b，c，d，e，f)T；

進(jìn)而求得二階梯度為

式中：x0、y0為單位方向矢量。

由于采用了與計(jì)算三角單元具有共用頂點(diǎn)的所有單元進(jìn)行計(jì)算，實(shí)際計(jì)算過程中式（10）一般為超定問題；如果式（10）欠定，可進(jìn)一步向周邊擴(kuò)展計(jì)算集合，以滿足最小二乘計(jì)算要求。

1.3 目標(biāo)函數(shù)最優(yōu)化求解

反演目標(biāo)函數(shù)的求解為最優(yōu)化問題，對于反演參數(shù)數(shù)量不大的一維問題可以直接采用奇異值分解；而對于二維或三維反問題常采用優(yōu)化解法，如非線性共軛梯度法、高斯-牛頓法、擬牛頓法等。由于高斯-牛頓求解效率高，本次研究采用高斯-牛頓法求解。將式（3）、（7）代入式（1），可得第k次高斯-牛頓迭代公式為

式中：

其中，JMT為AMT數(shù)據(jù)對模型參數(shù)的偏導(dǎo)數(shù)矩陣；JDC為DC數(shù)據(jù)對模型參數(shù)的偏導(dǎo)數(shù)矩陣。偏導(dǎo)數(shù)矩陣是滿陣，需要很大的存儲量，而且其計(jì)算往往也是反問題求解中最耗時(shí)的部分。為有效地減少靈敏度矩陣計(jì)算的工作量，通常采用互換定理計(jì)算偏導(dǎo)數(shù)矩陣；基于伴隨方程計(jì)算靈敏度矩陣的計(jì)算量與互換定理相當(dāng)，采用均勻半空間解析解作為伴隨方程解的近似方法可大大減少計(jì)算量。研究工作采用迭代法解式（12），應(yīng)用了一種隱式的偏導(dǎo)數(shù)矩陣求解方案，將偏導(dǎo)數(shù)矩陣及其轉(zhuǎn)置與向量的乘積轉(zhuǎn)換為正演計(jì)算，避免顯式存儲偏導(dǎo)數(shù)矩陣與海森矩陣，進(jìn)而有效地減少內(nèi)存需求。式（12）中，數(shù)據(jù)變化量

為第k次正演得到的計(jì)算數(shù)據(jù)與觀測數(shù)據(jù)的差。為確保每一次迭代目標(biāo)函數(shù)可以充分下降，對模型更新步長進(jìn)行線性搜索，有

χ為線性搜索步長，滿足

則

求解式（12），由于右端矩陣條件數(shù)大，直接解法很難求得穩(wěn)定解。采用正交分解方法進(jìn)行求解可提高反演穩(wěn)定性，但其計(jì)算量及內(nèi)存需求較大；采用共軛梯度法（CG）、廣義最小余量法（GMRES）等迭代解法進(jìn)行不精確求解，穩(wěn)定性較好；由CG發(fā)展而來的雙共軛梯度穩(wěn)定算法（BICGSTAB），計(jì)算效率高、穩(wěn)定性好，本文采用該算法求解高斯-牛頓方程。

1.4 漸進(jìn)反演算法

正則化技術(shù)依賴于先驗(yàn)信息，當(dāng)先驗(yàn)信息不足或者不正確時(shí)將會(huì)產(chǎn)生錯(cuò)誤的反演結(jié)果［21］，減少反演網(wǎng)格單元的數(shù)量與提高反演單元質(zhì)量是實(shí)現(xiàn)穩(wěn)定反演的重要途徑。生成高質(zhì)量網(wǎng)格的難度在于如何選擇單元細(xì)化或組合的標(biāo)準(zhǔn)，最直接的方法是通過分析以下模型分辨率矩陣［15］來實(shí)現(xiàn)，然而對于二維和三維反問題，由于其計(jì)算量過大，難于實(shí)現(xiàn)。

研究工作利用模型靈敏度信息，進(jìn)行自適應(yīng)網(wǎng)格生成，定義模型靈敏度為

模型靈敏度Sj是模型mj改變時(shí)對整個(gè)觀測數(shù)據(jù)集的影響。不考慮數(shù)據(jù)方差與穩(wěn)定因子的條件下，Sj實(shí)質(zhì)上反映了高斯-牛頓優(yōu)化過程中海森矩陣的主對角元素。通過Sj對網(wǎng)格進(jìn)行優(yōu)化實(shí)質(zhì)是調(diào)整海森矩陣的主對角元素，從而改善海森矩陣性質(zhì)，進(jìn)而得到高質(zhì)量的反演網(wǎng)格。

基于研發(fā)的網(wǎng)格優(yōu)化算法，進(jìn)一步研究了一種變網(wǎng)格反演技術(shù)。變網(wǎng)格算法一般采用兩種思路，一是先產(chǎn)生粗網(wǎng)格然后通過物性梯度［22］、多尺度算法［23］逐步細(xì)化；另一種則是先產(chǎn)生細(xì)網(wǎng)格，然后組合具有相似物性的單元［24］。此外，還有學(xué)者研究了基于角點(diǎn)與邊檢測技術(shù)的智能網(wǎng)格反演方法［25］。本文采用的漸進(jìn)網(wǎng)格反演以模型靈敏度為依據(jù)進(jìn)行網(wǎng)格細(xì)化，反演由粗網(wǎng)格到細(xì)網(wǎng)格逐步進(jìn)行，反演流程如圖2所示。計(jì)算過程如下：

Step 1讀入數(shù)據(jù)，設(shè)置漸進(jìn)網(wǎng)格反演迭代最大次數(shù)、網(wǎng)格優(yōu)化比例、最小單元面積等信息。

Step 2從數(shù)據(jù)信息生成正演模型的幾何描述，設(shè)定反演區(qū)域。

Step 3進(jìn)行非結(jié)構(gòu)化三角單元剖分。

Step 4當(dāng)前網(wǎng)格開始最優(yōu)化反演。設(shè)定當(dāng)前網(wǎng)格高斯-牛頓反演迭代次數(shù)、BICGSTAB迭代次數(shù)、收斂條件、正則化因子初值，設(shè)置初始模型和參考模型（初始網(wǎng)格采用均勻半空間值，后面采用前一重網(wǎng)格反演結(jié)果作為初始與參考模型）。

Step 5高斯-牛頓反演迭代，求取模型改變量。

Step 6線性搜索，并更新模型，調(diào)整正則化因子。

Step 7當(dāng)前網(wǎng)格進(jìn)行高斯-牛頓反演終止條件判斷，滿足進(jìn)入下一步，不滿足返回Step 5。

Step 8漸進(jìn)網(wǎng)格反演終止條件判斷，滿足退出反演，不滿足進(jìn)入Step 9。

Step 9互換定理計(jì)算雅克比矩陣，并計(jì)算模型靈敏度Sj。

Step 10分析模型靈敏度的平均值、方差，選擇待優(yōu)化單元。

Step 11返回Step 3，生成優(yōu)化網(wǎng)格。

圖2 漸進(jìn)網(wǎng)格反演流程Fig.2 Flowchart of step by step inversion

漸進(jìn)反演從粗網(wǎng)格開始，通過減少單元數(shù)量來減少反演對正則化因子的依賴；網(wǎng)格細(xì)化后，以前一重網(wǎng)格反演得到的模型作為下次網(wǎng)格的初始模型與參考模型，讓反演逐步進(jìn)行，進(jìn)一步保障反演的穩(wěn)定性，同時(shí)避免反演過程陷入局部極小或其他偽解。網(wǎng)格細(xì)化以模型靈敏度為依據(jù)進(jìn)行局部細(xì)化，具體算法：對模型靈敏度進(jìn)行統(tǒng)計(jì)分析，按設(shè)定比例選擇模型靈敏度大的單元，本文研究取反演單元數(shù)量的10%進(jìn)行優(yōu)化；同時(shí)，為了避免過度細(xì)化網(wǎng)格，通過設(shè)定最小單元面積進(jìn)行約束，即只對三角單元面積大于預(yù)設(shè)值的單元進(jìn)行細(xì)化。

2 理論模型算例

設(shè)計(jì)如圖3所示的二維模型，背景電阻率為100 Ω·m，共設(shè)置38個(gè)異常體，異常體均為截面為矩形的四棱柱，低阻體電阻率為10 Ω·m，高阻體電阻率為1 000 Ω·m。其中，接近地表交替有14個(gè)低阻和14個(gè)高阻異常體，異常體截面寬10 m、高30 m、上頂埋深為20 m，相鄰兩個(gè)異常體間隔為100 m；稍深交替有4個(gè)低阻與3個(gè)高阻體，異常體截面寬100 m、高100 m、上頂埋深為100 m，相鄰兩個(gè)異常體間隔為400 m；再向下，位于斷面中部有1個(gè)低阻體，異常體截面寬400 m、高200 m、上頂埋深為300 m；最深部有1個(gè)低阻和1個(gè)高阻異常體，異常體截面寬500 m、高500 m，上頂埋深為500 m，兩個(gè)異常體間隔為600 m。

近地表異常體用于模擬AMT的靜態(tài)效應(yīng)，靜態(tài)效應(yīng)嚴(yán)重干擾視電阻率，通常采用校正方法進(jìn)行壓制［26-27］；為避免校正方法引起的數(shù)據(jù)誤差，可通過其他淺層方法進(jìn)行補(bǔ)充［28-31］，進(jìn)而開展綜合解釋和聯(lián)合反演以提高解釋精度。

圖3 理論模型Fig.3 Synthetic model

直流電阻率法采用二極裝置，共布設(shè)電極320根，相鄰電極距10 m，最大極距200 m，數(shù)據(jù)集點(diǎn)數(shù)6 000個(gè)；AMT測點(diǎn)距40 m，共80個(gè)測點(diǎn)，測量頻段1～4 096 Hz，按2的指數(shù)共13個(gè)頻點(diǎn)，將正演數(shù)據(jù)加入3%隨機(jī)噪聲作為反演數(shù)據(jù)集。

圖4為正演得到的TM模式卡尼亞電阻率與相位。二維條件下TM模式卡尼亞電阻率受靜態(tài)效應(yīng)影響嚴(yán)重，斷面圖出現(xiàn)條帶狀異常，無法識別深部異常體；靜態(tài)效應(yīng)對相位影響小，相較視電阻率斷面圖可以分辨更多的異常體。

采用本文介紹的漸進(jìn)反演算法，分別對DC、AMT、DC與AMT數(shù)據(jù)進(jìn)行漸進(jìn)反演。反演過程對網(wǎng)格進(jìn)行5次細(xì)化，反演結(jié)果如圖5所示。

圖4 理論模型TM模式正演結(jié)果Fig.4 TM-mode synthetic data

圖5 理論模型反演斷面Fig.5 Sections of synthetic data inversion

圖5 a為直接反演DC數(shù)據(jù)所得斷面，反演結(jié)果沒有得到深部異常體；近地表的28個(gè)異常體反演效果較好，但受稍深的7個(gè)異常體影響存在畸變。圖5b為AMT反演結(jié)果，從反演斷面可以判斷出幾乎所有異常體，但是很多異常相互聯(lián)結(jié)無法分離；近地表異常反演效果不理想，異常體位置有高、低阻體出現(xiàn)，但反演所得電阻率值與真實(shí)值相差較大，說明模擬設(shè)定的點(diǎn)距與頻率條件下AMT對淺層不均勻體有一定的反演能力。圖5c為聯(lián)合反演兩個(gè)數(shù)據(jù)集的反演結(jié)果，反演有效地給出了淺層與深部所有異常的位置；相比于單獨(dú)反演，近地表28個(gè)異常體的電阻率與空間位置均較單獨(dú)反演更接近真實(shí)值；稍深的7個(gè)異常體形態(tài)與設(shè)定模型有一定的出入，但異常體的位置與反演電阻率值均較單獨(dú)反演精確；只有聯(lián)合反演可分辨斷面中部獨(dú)立的低阻異常體；深部低阻異常的反演效果優(yōu)于高阻，分析原因?yàn)樵诒疚乃捎玫哪Ｐ驮O(shè)置下，低阻二度體相較高阻二度體更易形成沿走向方向的電流場，有利于反演。

圖6 反演區(qū)網(wǎng)格變化Fig.6 Mesh of inversion domain

圖7 正則化因子、數(shù)據(jù)均方根誤差、模型誤差變化曲線Fig.7 Curves of optimum multiplier,RMS,and stabilizer

表1 漸進(jìn)反演反演次數(shù)、RMS與模型誤差統(tǒng)計(jì)表Tab.1 Iteration times,RMS,and model error of AMT inversion

為進(jìn)一步分析漸進(jìn)反演過程，現(xiàn)對AMT反演進(jìn)行分析，反演過程網(wǎng)格變化如圖6所示。表1給出了各次反演網(wǎng)格單元數(shù)量、總單元數(shù)量，反演迭代開始和結(jié)束時(shí)的數(shù)據(jù)均方根誤差、模型誤差值。圖7為正則化因子、數(shù)據(jù)均方根誤差（RMS）、模型誤差變化曲線。反演區(qū)單元數(shù)從836到1 932，逐漸增加，第1、2次網(wǎng)格反演，RMS值快速下降，后面優(yōu)化網(wǎng)格RMS值下降量逐漸減少，趨于穩(wěn)定；反演過程中，初始網(wǎng)格正則化因子出現(xiàn)增大現(xiàn)象，隨后穩(wěn)定下降，表明漸進(jìn)網(wǎng)格方法對正則化因子的初始依賴不大，可以通過自動(dòng)調(diào)整確保反演穩(wěn)定進(jìn)行；最后兩重網(wǎng)格數(shù)據(jù)誤差下降不多，更多的是對模型進(jìn)行優(yōu)化；由于地表的靈敏度較大，網(wǎng)格總體由上到下逐漸加密，從圖6c～6f看，由電阻率分布引起的網(wǎng)格變化過程清晰，體現(xiàn)了基于模型靈敏度的網(wǎng)格優(yōu)化特點(diǎn)。

3 實(shí)測數(shù)據(jù)反演

采用本文的反演算法，對某沉積巖區(qū)的DC與AMT勘探數(shù)據(jù)開展聯(lián)合反演研究，并與單獨(dú)反演DC與AMT的結(jié)果進(jìn)行對比分析?？碧狡拭骈L度3 500 m，其中直流電阻率法采用10 m電極距、溫納α裝置測量，最大供電極距900 m，共采集有效直流電阻率測點(diǎn)7 474個(gè)；AMT采用50 m點(diǎn)距，頻率范圍1～9 600 Hz，共41個(gè)頻點(diǎn)，共采集有效測深點(diǎn)70個(gè)。

反演結(jié)果如圖8所示，反演斷面電阻率變化范圍1～1 000 Ω·m，主體電阻率不高。DC反演斷面，測線前600 m表現(xiàn)為中間高的三層結(jié)構(gòu)，且近地表存在不連續(xù)低阻；而600 m以后，出現(xiàn)連續(xù)層狀低阻體，其上部存在連續(xù)性較差的高阻，近地表出現(xiàn)不連續(xù)低阻；考慮到所采用的工作裝置最大供電極距為900 m，勘探深度有限，300 m以下與初始模型電阻率相差不大的區(qū)域應(yīng)當(dāng)達(dá)到了直流反演的最大深度。AMT反演斷面，在測線起點(diǎn)端左下方出現(xiàn)了一個(gè)中高阻體，并向淺部發(fā)展，在接近水平坐標(biāo)500 m的位置接近地表；近地表500～680 m之間出現(xiàn)近直立低阻與高阻相伴異常；680 m后，電阻率性質(zhì)與DC反演斷面分布形態(tài)相似，但AMT反演斷面，低阻的分布區(qū)明顯比直流更大，這與直流的勘探深度不足有關(guān)。相較DC與AMT單獨(dú)反演斷面，聯(lián)合反演得到了更為合理的底部高阻，高阻頂界面相比AMT變化緩，與實(shí)質(zhì)地質(zhì)條件相符。100～600 m深度出現(xiàn)了更精細(xì)的結(jié)構(gòu)，分析原因?yàn)椋孩佼?dāng)極距增大時(shí)直流溫納α裝置本身更適合反演層狀結(jié)構(gòu)。②工區(qū)淺部電阻率低，電磁法的趨膚深度小、垂向分辨率高。③AMT本身的橫向分辨率較高，所以聯(lián)合反演取得了理想的淺層識別能力；水平坐標(biāo)500～680 m近直立低阻與高阻相伴異常形態(tài)更清晰；測線起點(diǎn)淺部高阻體邊界更為明顯。

圖8 實(shí)測數(shù)據(jù)反演斷面Fig.8 Sections of field data inversion

4 結(jié)論

為提高反演精度，本文開展了DC和AMT自適應(yīng)漸進(jìn)聯(lián)合反演研究。通過對反演方法的研究和模型試算，取得以下幾點(diǎn)認(rèn)識：

（1）非結(jié)構(gòu)化三角網(wǎng)格剖分適用于直流電阻率與音頻大地電磁兩種空間采樣與分辨率差別較大的方法進(jìn)行聯(lián)合反演時(shí)對剖分的需求。

（2）基于模型靈敏度的漸進(jìn)網(wǎng)格反演算法，減小了正則化方法對穩(wěn)定因子的依賴，提高了反演的穩(wěn)定性，同時(shí)減少了對正則化因子進(jìn)行搜索的計(jì)算量。

（3）聯(lián)合反演相較單一方法提高了反演模型的分辨率，可實(shí)現(xiàn)AMT產(chǎn)生靜態(tài)效應(yīng)異常體的直接反演。