二維直流電阻率與音頻大地電磁自適應漸進聯合反演

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

（1.同濟大學海洋與地球科學學院，上海200092；2.東華理工大學地球物理與測控技術學院，江西南昌330013；3.中國地質調查局南京地質調查中心，江蘇南京210016）

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

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

1 DC與AMT正則化聯合反演

為更好地擬合地表與地下地質單元，實現野外工作點距差別大、勘探深度與分辨率不同的AMT與DC兩種方法的高效、高精度正演，本文采用非結構化三角網格進行區域剖分；應用有限單元法進行DC與AMT正演［11-14］；采用基于最小結構模型的正則化方法構造非結構化三角網格條件下的模型二階平滑計算方法，采用高斯-牛頓最優化方法求解反演目標函數，并通過雙共軛梯度穩定算法（BICGSTAB）不完全求解高斯-牛頓方程，以確保反演過程穩定，同時減少反演對內存的需求。

1.1 正則化目標函數

地球物理反問題通常為病態問題，正則化是穩定求解反問題的有效手段［15-16］。正則化反演的目標函數為

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

DC與AMT聯合反演的變量為介質電阻率，采用以下變換函數［17-18］以確保反演電阻率在實際巖、礦石的取值范圍內：

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

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

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

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

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

式中：βs，βx，βy為模型誤差三部分之間的比例系數；mapr為先驗模型；模型誤差可統一表示為

式中：Wm為統一定義形式的模型誤差矩陣。

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

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

1.2 非結構化三角網格模型梯度計算

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

圖1 非結構三角網格梯度Fig.1 Gradient calculation of unstructured triangular mesh

假設物性分布為空間坐標的二次函數，有

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

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

進而求得二階梯度為

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

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

1.3 目標函數最優化求解

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

式中：

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

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

χ為線性搜索步長，滿足

則

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

1.4 漸進反演算法

正則化技術依賴于先驗信息，當先驗信息不足或者不正確時將會產生錯誤的反演結果［21］，減少反演網格單元的數量與提高反演單元質量是實現穩定反演的重要途徑。生成高質量網格的難度在于如何選擇單元細化或組合的標準，最直接的方法是通過分析以下模型分辨率矩陣［15］來實現，然而對于二維和三維反問題，由于其計算量過大，難于實現。

研究工作利用模型靈敏度信息，進行自適應網格生成，定義模型靈敏度為

模型靈敏度Sj是模型mj改變時對整個觀測數據集的影響。不考慮數據方差與穩定因子的條件下，Sj實質上反映了高斯-牛頓優化過程中海森矩陣的主對角元素。通過Sj對網格進行優化實質是調整海森矩陣的主對角元素，從而改善海森矩陣性質，進而得到高質量的反演網格。

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

Step 1讀入數據，設置漸進網格反演迭代最大次數、網格優化比例、最小單元面積等信息。

Step 2從數據信息生成正演模型的幾何描述，設定反演區域。

Step 3進行非結構化三角單元剖分。

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

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

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

Step 7當前網格進行高斯-牛頓反演終止條件判斷，滿足進入下一步，不滿足返回Step 5。

Step 8漸進網格反演終止條件判斷，滿足退出反演，不滿足進入Step 9。

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

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

Step 11返回Step 3，生成優化網格。

圖2 漸進網格反演流程Fig.2 Flowchart of step by step inversion

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

2 理論模型算例

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

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

圖3 理論模型Fig.3 Synthetic model

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

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

采用本文介紹的漸進反演算法，分別對DC、AMT、DC與AMT數據進行漸進反演。反演過程對網格進行5次細化，反演結果如圖5所示。

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

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

圖5 a為直接反演DC數據所得斷面，反演結果沒有得到深部異常體；近地表的28個異常體反演效果較好，但受稍深的7個異常體影響存在畸變。圖5b為AMT反演結果，從反演斷面可以判斷出幾乎所有異常體，但是很多異常相互聯結無法分離；近地表異常反演效果不理想，異常體位置有高、低阻體出現，但反演所得電阻率值與真實值相差較大，說明模擬設定的點距與頻率條件下AMT對淺層不均勻體有一定的反演能力。圖5c為聯合反演兩個數據集的反演結果，反演有效地給出了淺層與深部所有異常的位置；相比于單獨反演，近地表28個異常體的電阻率與空間位置均較單獨反演更接近真實值；稍深的7個異常體形態與設定模型有一定的出入，但異常體的位置與反演電阻率值均較單獨反演精確；只有聯合反演可分辨斷面中部獨立的低阻異常體；深部低阻異常的反演效果優于高阻，分析原因為在本文所采用的模型設置下，低阻二度體相較高阻二度體更易形成沿走向方向的電流場，有利于反演。

圖6 反演區網格變化Fig.6 Mesh of inversion domain

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

表1 漸進反演反演次數、RMS與模型誤差統計表Tab.1 Iteration times,RMS,and model error of AMT inversion

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

3 實測數據反演

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

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

圖8 實測數據反演斷面Fig.8 Sections of field data inversion

4 結論

為提高反演精度，本文開展了DC和AMT自適應漸進聯合反演研究。通過對反演方法的研究和模型試算，取得以下幾點認識：

（1）非結構化三角網格剖分適用于直流電阻率與音頻大地電磁兩種空間采樣與分辨率差別較大的方法進行聯合反演時對剖分的需求。

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

（3）聯合反演相較單一方法提高了反演模型的分辨率，可實現AMT產生靜態效應異常體的直接反演。