優化策略的二維大地電磁光滑聚焦反演研究

白寧波 周君君 胡祥云*

(①中國地質大學(武漢)地質探測與評估教育部重點實驗室，湖北武漢 430074；②中國地質大學(武漢)地球物理與空間信息學院，湖北武漢 430074)

0 引言

大地電磁測深(magnetotelluric，MT)作為發展較早的地球物理勘探方法之一，廣泛應用于地球固體礦產和石油天然氣的勘探等領域。當前大地電磁反演已經發展到三維成像階段，然而三維反演的計算量較大，對硬件要求高，因此二維反演在實際應用中仍然具有很大優勢[1]。典型的二維反演有OCCAM反演[2-3]、快速松弛反演(RRI)[4]、簡化基OCCAM反演(REBOCC)[5]及非線性共軛梯度反演(NLGG)[6]等。盡管這些反演方法都可以得到很好的反演效果，但分別是以最小模型約束、最平緩模型約束或最光滑模型約束作為穩定泛函，因此只能得到平滑的反演效果，而不能得到清晰的地質體分界面。針對這個問題，Last等[7]首先提出用最小支撐泛函(minimum support，MS)作為穩定泛函，提高了塊狀結構的分辨率。隨后，Portniaguine等[8]在此基礎上提出使用最小梯度支撐(minimum gradient support，MGS)作為穩定泛函進行地球物理反演，得到了清晰的地質體分界面。Zhang等[9]和Zhdanov[10]針對MGS做了進一步研究，取得了較好的聚焦反演效果；張羅磊等[11]結合MGS穩定泛函和OCCAM方法，反演結果突出了對尖銳電性邊界的刻畫。隨著反演理論的發展，一些新的穩定泛函被引入，如Sun等[12]提出的反正切穩定泛函、Hu等[13]提出的反余切穩定泛函、Zhao等[14]和Hu等[15]提出的指數型穩定泛函，以及Xiang等[16]提出的最小支撐梯度穩定泛函(minimum support gradient，MSG)。盡管基于這些穩定泛函可以得到清晰的地質體分界面，但是聚焦反演可能會使構造形態發生畸變，反演結果不準確。

采用高斯—牛頓法求解反演目標泛函時，由于正則化因子是數據擬合泛函和模型穩定泛函的折中參數，正則化因子過大會過于強調模型穩定泛函，導致數據總體欠擬合；反之，則容易產生虛假的反演構造。因此，正則化因子的選取對反演結果影響很大。關于正則化因子選取的典型方法，主要有Hansen等[17]提出的L曲線法、吳小平等[18]提出的自適應遞減方法及陳小斌等[19]提出的完全自適應正則化方法(CMD方案)。此外，也有一些改進的正則化因子方案，如張羅磊等[11]、樸英哲等[20]及向陽等[21]都對正則化因子的選取提出了改進策略，這在一定程度上提高了反演效率。

本文對光滑反演和聚焦反演的缺陷進行分析，進而提出一種新的反演目標泛函。該目標泛函將最光滑約束與MSG穩定泛函進行加權結合，利用高斯—牛頓法對目標泛函進行求解，不僅可以得到穩定的反演結果，還可以使地質體分界面變得更加清晰。同時，在進行反演迭代時，本文提出利用Nelder-Mead優化算法優化Morozov偏差原理選取正則化因子的方法，不僅彌補了Morozov 偏差原理后期收斂速度慢的不足，還加快了反演算法的收斂速度。同時，本文采用自適應衰減的聚焦因子。文中對典型模型進行反演，通過對比不同反演策略的反演結果，驗證了本文算法的優勢。最后，針對山西陽高縣的一條實測數據進行反演，結果驗證了光滑聚焦反演的有效性和可靠性。

1 理論分析

1.1 反演理論

根據正則化理論，大地電磁反演目標泛函表達式為

P(m,d)=f(m,d)+αs(m)

(1)

式中：P(m,d)是參數目標泛函，m=[m1,m2,…,mM]T為模型參數向量，其中M是模型參數的數量，d=[d1,d2,…,dN]T是觀測數據向量，N是觀測數據的數量;α是正則化參數；f(m,d)和s(m)分別是數據目標泛函和模型目標泛函。f(m,d)的表達式為

(2)

式中：Wd是數據加權矩陣；F是正演算子。

采用不同的模型目標泛函進行約束，會得到不同的反演結果。目前常用的三種穩定泛函是最小模型約束、最平緩模型約束和最光滑模型約束，這三種約束泛函盡管可以得到較好的反演結果，卻不能得到清晰的地質體分界面。聚焦反演的引入很好地解決了這一難題。另一方面，采用聚焦反演很可能導致反演結果嚴重聚焦。為了避免這一缺陷，本文采用最光滑約束和MSG同時對模型目標泛函進行約束，具體表達式為

s(m)=γs1(m)+(1-γ)s2(m)

(3)

(4)

(5)

式(4)中的MSG穩定泛函s1(m)是在Last等[7]提出的MS泛函的基礎上進一步求空間梯度得到的，Zhdanov[22]已證明MS滿足正則穩定器的Tikhonov準則，可以作為穩定泛函對模型進行約束。因此，MSG同樣也可以作為模型穩定泛函，Zhao等[14]對此作了相關描述。

式(4)中，聚焦參數β過大會導致聚焦效果不明顯，過小則可能使泛函產生奇異，所以選擇合適的聚焦因子對反演結果同樣重要。針對β的選取，Zhdanov等[23]提出了L曲線法，但每次迭代都要進行曲率計算。為此，本文提出下列自適應衰減因子

β=e-ωk

(6)

式中：ω為控制系數，本文取1；k為當前的迭代次數，k=1,2,…，K，K表示停止迭代次數。

因此式(3)和式(4)可寫成泛函形式

(7)

(8)

式中We=WW，W的表達式為

則式(1)中的參數目標泛函可以轉化為

(9)

對于模型加權矩陣We，本文采用Portniaguine等[8]提出的方法進行處理。假設第k次迭代時，當前模型參數為mk,則令We=We(mk)，這樣在每次迭代過程中，模型加權矩陣We就可當做一個常數矩陣。

對式(9)進行一階泰勒展開，并令目標泛函的一階變分等于零。為保證迭代的穩定性，對反演迭代公式先取對數后再反演，則可得到高斯—牛頓法的更新迭代公式

lnmk+1-lnmapr={(WdJk)HWdJk+

(10)

Jk(lnmk-lnmapr)]

(11)

式中：J為雅克比矩陣；H表示共軛轉置。

由式(10)可知，正則化因子α對反演是否收斂至關重要。下文將闡述如何基于Nelder-Mead優化算法的Morozov偏差原理確定正則化因子。

1.2 Morozov偏差原理

利用Morozov偏差原理確定正則因子的應用詳見文獻[24]。該原理決定了正則因子存在反演后期收斂速度變慢甚至不收斂的缺陷。Morozov偏差原理是根據后驗策略選擇最優的正則化參數，而且可以自動選擇唯一的正則化參數。假設數據采集時誤差水平已知，則需要滿足以下偏差方程

(12)

φ(α)≈φ(αk)+φ'(αk)(α-αk)

(13)

令φ(α)=0，可以得到牛頓法的迭代格式

(14)

采用Morozov偏差原理決定了正則化因子盡管能使反演結果收斂，但是算法后期的收斂速度很慢，反演效率低，因此需要對Morozov偏差原理的結果做進一步優化處理。

1.3 Nelder-Mead優化算法

Nelder-Mead優化算法又稱為可變的多面體搜索法，是一種無約束的直接搜索方法。該算法的基本思想是首先利用起始點x0構建一個具有n+1個頂點的線性多面體(x0,x1,…，xn)，通過對比各個頂點的目標函數值確定各頂點的優劣。假設計算可得xL為最差點，xH為最優點，xC為次優點，可通過啟發式的反射、擴張、壓縮等運算得到新的頂點(xl，xm，xn)，然后用較好的頂點替換最差的頂點xL組成新的多面體。Nelder-Mead優化算法如圖1所示，Yildiz等[25]對Nelder-Mead優化算法的具體過程進行了描述。最后，對該流程反復進行迭代運算，最終逼近目標函數的最優解。Nelder-Mead優化算法收斂速度快，對局部的搜索能力強，可用極小的時間成本搜尋局部最優值。

圖1 Nelder-Mead優化算法示意圖

反演的均方根誤差RMS的計算公式為

(15)

本文基于優化策略的光滑聚焦反演的具體步驟如下。

(1)給出初始模型m0、先驗模型mapr、初始正則化因子α1、權重因子γ、反演迭代次數K、目標擬合差及采用Nelder-Mead優化算法尋優的最大迭代次數N。

(2)計算數據的RMS，判斷是否滿足給定的停止條件。若滿足，則停止迭代；若不滿足，則轉步驟(3)。

(3)計算當前模型的雅克比矩陣Jk、F(mk)、W和We(mk)。利用Nelder-Mead優化算法對正則化因子αk進一步優化，得到新的正則化因子α'k。利用式(10)更新模型。

(4)采用Morozov偏差原理決定下一次迭代的正則化因子αk+1，然后令k=k+1，轉至步驟(2)。

1.4 矢量化與并行

進行大地電磁二維反演時需要反復調用正演子程序和雅可比矩陣的計算子程序，考慮到反演效率，反演時本文采用矢量化思想和Matlab的并行策略進行編程。矢量化編程是利用單元網格局部編碼和整體編碼的策略，將整體的剛度矩陣中的非零行和非零列直接寫成了矩陣形式，通過一次性賦值即可得到整體的剛度矩陣，這樣可避免使用多次循環進行賦值，提高了賦值效率，減少了正演的時間。Matlab并行計算是利用正演時各個頻率間的正演過程相互獨立，調用CPU多核進行并行計算，大大提高了正演速度。以上兩種策略可保證較快的反演過程。

2 模型試驗

2.1 Sasaki模型

為了驗證算法的優勢和可靠性，選用典型的Sasaki模型[26]進行反演，見圖1。背景模型是電阻率為50Ω·m的均勻半空間，包含不同的高阻和低阻異常體。測點位于地面，間距為0.3km，沿y軸分布，分布區域為-9～9km，頻點數為25，頻率范圍是0.1～100Hz，按對數等間距分布。

反演過程中，對于初始模型m0選取50Ω·m的均勻半空間，設mapr=m0，目標擬合差設定為0.01，反演迭代次數為10，對正則化因子進行優化時，Nelder-Mead優化算法的迭代次數N取10。

為了對比分析，采用遞減策略選取正則化因子的OCCAM進行反演，結果見圖3a；圖3b為優化策略的OCCAM反演結果，即式(3)中γ=0時的反演結果；圖3c～圖3f為γ分別取0.1、0.5、0.9和1.0的反演結果，其中γ=1.0即對應MSG反演。圖4為對應圖3的不同方法擬合差迭代曲線。由圖3和圖4可見，優化策略的反演收斂速度更快，所需時間更少。對比圖3a、圖3b與圖3c～圖3f可知，目標函數含有MSG項的穩定泛函可反演得到更清晰的地質體界面。從圖3c～圖3f還可見，隨著γ的增加，地質體的分界面越來越清晰，聚焦效果也越來越明顯，但低阻體電阻率偏離真值的程度也越來越高。圖5是對應圖3a～圖3f的反演結果與真實模型的絕對誤差剖面。對比圖3與圖5可知，γ=0.9時的反演效果相對最好。

圖2 Sasaki 模型示意圖

圖3 不同策略的OCCAM電阻率反演結果

圖4 對應圖3的不同策略反演收斂曲線對比

圖5 對應圖3的不同反演策略反演結果與真實值的絕對誤差

進行光滑聚焦聯合反演時，模型泛函s1(m)中聚焦因子β的選取對聚焦效果有很大影響。因此，計算γ=0.9時，β分別取1.0、0.5和0.1的反演結果見圖6。由圖可見，隨著β值的降低，聚焦效果越來越明顯，因此進行聚焦反演時需分析選取合適的聚焦因子。而式(6)因采用了自適應遞減的聚焦因子，故不必考慮這一問題。

圖6 優化策略聚焦反演時不同β值的反演結果對比(γ=0.9)

2.2 楔形模型

構建一個二維楔形模型，驗證本文方法反演結果的可靠性和精度性，模型參數見圖7。設置測點間距為0.2km。圖8是γ=0、0.9時的反演結果及其與模型真實電阻率的誤差分布，圖9是γ=0.9時的迭代收斂曲線。從圖8電阻率反演結果可以看出，光滑聚焦反演結果中地質體界面更加清晰；從圖9收斂曲線可以看出，對于楔形模型，采用Nelder-Mead優化算法的Morozov偏差原理確定正則化因子仍有很好的收斂效果。

圖7 二維楔形模型示意圖

圖8 二維楔形模型優化策略聚焦反演結果(上)及其與模型真實電阻率的差值(下)

圖9 二維楔形模型優化策略(γ=0.9)聚焦反演收斂曲線

2.3 實測數據反演

為了進一步檢驗算法的有效性和可靠性，對山西省陽高縣某地熱勘探測線的實測數據進行了反演。該測線位于陽高縣城北側大同盆地北緣，測線經過一個斷裂，該斷裂是在燕山期“陽高破碎帶”的基礎上繼承形成的。測點點距為500m，共22個測點，所測數據為MT。圖10a是該測線實測數據的擬斷面圖，圖10b是采用有限差分正演和非線性共軛梯度(商業軟件Winglink)的反演結果，圖10c是本文光滑聚焦的電阻率反演結果。從圖10b和圖10c可以看出，在深度2～6km范圍內都有一個明顯的低阻體，根據已有的實際資料推斷，低阻區域為斷裂含水帶，是斷裂構造作用造成巖石破碎，導致透水性增大形成的。對比圖10b與圖10c可見，光滑聚焦反演的結果與Winglink軟件反演的結果大致相同，驗證了本文算法的有效性和可靠性。由圖10c還可以看到，光滑聚焦反演的低阻區域較為明顯，并且對低阻體邊界的刻畫較Winglink軟件反演結果更清晰，因此光滑聚焦反演算法對于實際大地電磁資料的反演也是有效的。

圖10 實測數據及電阻率反演結果

3 結論

(1)本文提出了一種新的目標泛函，利用高斯—牛頓法求解反演目標函數，實現了二維大地電磁數據的穩定聚焦反演，可以得到清晰的地質界面。

(2)針對反演迭代過程，提出了利用Nelder-Mead優化算法優化Morozov偏差原理確定正則化因子的優化策略，減少了反演迭代次數，加快了反演收斂速度。

(3)通過對典型模型和實測數據的反演分析，驗證了光滑聚焦反演的可靠性和有效性，為三維MT反演奠定了應用基礎。