白寧波 周君君 胡祥云*

(①中國(guó)地質(zhì)大學(xué)(武漢)地質(zhì)探測(cè)與評(píng)估教育部重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430074；②中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，湖北武漢 430074)

0 引言

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

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

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

1 理論分析

1.1 反演理論

根據(jù)正則化理論，大地電磁反演目標(biāo)泛函表達(dá)式為

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

(1)

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

(2)

式中：Wd是數(shù)據(jù)加權(quán)矩陣；F是正演算子。

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

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

(3)

(4)

(5)

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

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

β=e-ωk

(6)

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

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

(7)

(8)

式中We=WW，W的表達(dá)式為

則式(1)中的參數(shù)目標(biāo)泛函可以轉(zhuǎn)化為

(9)

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

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

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

(10)

Jk(lnmk-lnmapr)]

(11)

式中：J為雅克比矩陣；H表示共軛轉(zhuǎn)置。

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

1.2 Morozov偏差原理

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

(12)

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

(13)

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

(14)

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

1.3 Nelder-Mead優(yōu)化算法

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

圖1 Nelder-Mead優(yōu)化算法示意圖

反演的均方根誤差RMS的計(jì)算公式為

(15)

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

(1)給出初始模型m0、先驗(yàn)?zāi)Ｐ蚼apr、初始正則化因子α1、權(quán)重因子γ、反演迭代次數(shù)K、目標(biāo)擬合差及采用Nelder-Mead優(yōu)化算法尋優(yōu)的最大迭代次數(shù)N。

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

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

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

1.4 矢量化與并行

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

2 模型試驗(yàn)

2.1 Sasaki模型

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

反演過程中，對(duì)于初始模型m0選取50Ω·m的均勻半空間，設(shè)mapr=m0，目標(biāo)擬合差設(shè)定為0.01，反演迭代次數(shù)為10，對(duì)正則化因子進(jìn)行優(yōu)化時(shí)，Nelder-Mead優(yōu)化算法的迭代次數(shù)N取10。

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

圖2 Sasaki 模型示意圖

圖3 不同策略的OCCAM電阻率反演結(jié)果

圖4 對(duì)應(yīng)圖3的不同策略反演收斂曲線對(duì)比

圖5 對(duì)應(yīng)圖3的不同反演策略反演結(jié)果與真實(shí)值的絕對(duì)誤差

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

圖6 優(yōu)化策略聚焦反演時(shí)不同β值的反演結(jié)果對(duì)比(γ=0.9)

2.2 楔形模型

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

圖7 二維楔形模型示意圖

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

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

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

為了進(jìn)一步檢驗(yàn)算法的有效性和可靠性，對(duì)山西省陽(yáng)高縣某地?zé)峥碧綔y(cè)線的實(shí)測(cè)數(shù)據(jù)進(jìn)行了反演。該測(cè)線位于陽(yáng)高縣城北側(cè)大同盆地北緣，測(cè)線經(jīng)過一個(gè)斷裂，該斷裂是在燕山期“陽(yáng)高破碎帶”的基礎(chǔ)上繼承形成的。測(cè)點(diǎn)點(diǎn)距為500m，共22個(gè)測(cè)點(diǎn)，所測(cè)數(shù)據(jù)為MT。圖10a是該測(cè)線實(shí)測(cè)數(shù)據(jù)的擬斷面圖，圖10b是采用有限差分正演和非線性共軛梯度(商業(yè)軟件Winglink)的反演結(jié)果，圖10c是本文光滑聚焦的電阻率反演結(jié)果。從圖10b和圖10c可以看出，在深度2～6km范圍內(nèi)都有一個(gè)明顯的低阻體，根據(jù)已有的實(shí)際資料推斷，低阻區(qū)域?yàn)閿嗔押畮?，是斷裂?gòu)造作用造成巖石破碎，導(dǎo)致透水性增大形成的。對(duì)比圖10b與圖10c可見，光滑聚焦反演的結(jié)果與Winglink軟件反演的結(jié)果大致相同，驗(yàn)證了本文算法的有效性和可靠性。由圖10c還可以看到，光滑聚焦反演的低阻區(qū)域較為明顯，并且對(duì)低阻體邊界的刻畫較Winglink軟件反演結(jié)果更清晰，因此光滑聚焦反演算法對(duì)于實(shí)際大地電磁資料的反演也是有效的。

圖10 實(shí)測(cè)數(shù)據(jù)及電阻率反演結(jié)果

3 結(jié)論

(1)本文提出了一種新的目標(biāo)泛函，利用高斯—牛頓法求解反演目標(biāo)函數(shù)，實(shí)現(xiàn)了二維大地電磁數(shù)據(jù)的穩(wěn)定聚焦反演，可以得到清晰的地質(zhì)界面。

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

(3)通過對(duì)典型模型和實(shí)測(cè)數(shù)據(jù)的反演分析，驗(yàn)證了光滑聚焦反演的可靠性和有效性，為三維MT反演奠定了應(yīng)用基礎(chǔ)。