大地電磁二維DCG、NLCG和L-BFGS反演算法對(duì)比

李 磊，吳桂桔，談洪波

(中國(guó)地震局地震研究所 地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430071)

1 引 言

大地電磁法(Magnetotelluric，MT)作為地下深部結(jié)構(gòu)探測(cè)的一種重要手段，已經(jīng)被廣泛地應(yīng)用于地殼和上地幔的電性結(jié)構(gòu)探測(cè)工作[1,2]。目前，在所有地球物理方法中，MT是除地震方法以外探測(cè)殼幔結(jié)構(gòu)最重要的方法，這是因?yàn)樗粌H成本低，便于實(shí)際野外施工，而且分辨能力高，探測(cè)深度大。由于真實(shí)地下介質(zhì)的電性結(jié)構(gòu)常具有三維性特征，為了反演結(jié)果的真實(shí)性，MT野外實(shí)際采集資料最好采用三維反演方法[3,4]。雖然一些復(fù)雜的區(qū)域可能需要三維反演，但三維反演計(jì)算需要耗費(fèi)大量的時(shí)間，二維反演的計(jì)算時(shí)間更快(肖俊等，2020)[5]。一維反演可以作為二維或三維反演的初始模型(劉嘉文等，2018)[6]。此外，二維反演可以方便地加入地形的影響，這是一維反演所不具備的(岳明鑫等，2019)[7]。因此，MT二維反演仍是目前最主流的處理手段。

MT二維反演方法種類很多[8-13]，常用的方法包括Occam反演法[14-16]，快速松弛反演法(RRI)[17]， 非線性共軛梯度法(NLCG)[18]，數(shù)據(jù)空間Occam反演法(REBOCC)[19]以及數(shù)據(jù)空間共軛梯度法(DCG)[20]等。這些方法都屬于基于目標(biāo)函數(shù)梯度信息的最優(yōu)化迭代反演方法，但它們具有不同的特點(diǎn)，有的單次迭代計(jì)算速度快、內(nèi)存需求小，有的迭代收斂速度快。從理論上來(lái)看，在目標(biāo)函數(shù)中的模型約束項(xiàng)、數(shù)據(jù)擬合項(xiàng)采用統(tǒng)一定義的情況下，這些不同反演算法的結(jié)果應(yīng)該是基本一致的(在收斂條件也相同的情況下)。然而，實(shí)際情況卻并非如此，由于實(shí)際野外地質(zhì)結(jié)構(gòu)的復(fù)雜性[21]，實(shí)測(cè)MT資料使用不同的反演方法往往效果相差明顯。因此在實(shí)際應(yīng)用中，MT資料處理工作者都要面臨著如何選取更適合的反演方法的問(wèn)題。

本文將對(duì)DCG、NLCG以及有限內(nèi)存擬牛頓法(L-BFGS)[22]三種二維大地電磁反演方法進(jìn)行對(duì)比。雖然之前也有學(xué)者對(duì)MT二維反演算法進(jìn)行過(guò)對(duì)比[23]，但對(duì)比的主要是Occam和數(shù)據(jù)空間Occam類方法，沒(méi)有對(duì)DCG和L-BFGS對(duì)比過(guò)。本文首先介紹這幾種方法的基本理論，然后設(shè)計(jì)典型的地電模型，使用以上三種反演算法對(duì)這些模型進(jìn)行反演試算。本文使用TE、TM、Tzy三種不同模式下的數(shù)據(jù)進(jìn)行反演，同時(shí)對(duì)比不同反演算法的效果，以期得出不同反演方法對(duì)反演結(jié)果以及對(duì)地球物理解釋工作的影響，從而方便地球物理工作者在實(shí)際進(jìn)行大地電磁處理解釋工作時(shí)更好地選取MT二維反演方法。

2 大地電磁常用反演算法的基本理論

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

DCG、NLCG和 LBFGS三種MT二維反演方法都屬于正則化反演方法，它們的目標(biāo)函數(shù)可表示為:

φ=φd+αφm

(1)

式(1)中，φ是總目標(biāo)函數(shù);φm為模型約束目標(biāo)函數(shù)；φd是觀測(cè)數(shù)據(jù)目標(biāo)函數(shù)；α是正則化因子；m表示模型向量，d表述數(shù)據(jù)向量。

常用的模型約束函數(shù)有下列三種：

1)最小模型約束目標(biāo)函數(shù)，模型參數(shù)平方和達(dá)到最小值:

φm=‖m‖2=(m,m)min

(2)

2)最平緩模型約束目標(biāo)函數(shù)，模型參數(shù)的一階導(dǎo)的平方和達(dá)到最小值[2]:

φm=‖?m‖2=(?m,?m)min

(3)

3)最光滑模型約束目標(biāo)函數(shù)，模型參數(shù)的二階導(dǎo)的平方和達(dá)到最小值:

φm=‖?2m‖2=(?2m,?2m)min

(4)

2.2 數(shù)據(jù)空間的共軛梯度法(DCG)

DCG就是用共軛梯度法(CG)在數(shù)據(jù)空間來(lái)求解反演法方程，被稱為數(shù)據(jù)空間共軛梯度法[4]，所以先給出共軛梯度法的基本公式。

2.2.1 共軛梯度法(CG)

CG迭代過(guò)程：

(5)

(6)

搜索方向可以表達(dá)為：

(7)

(8)

(9)

Cd是數(shù)據(jù)協(xié)方差矩陣，Cm是模型協(xié)方差矩陣。所以梯度向量迭代方式可以改寫(xiě)成：

(10)

2.2.2 基于數(shù)據(jù)空間的共軛梯度法(DCG)

將上述CG算法引入到數(shù)據(jù)空間反演方法中，經(jīng)過(guò)推導(dǎo)可得到基于數(shù)據(jù)空間的共軛梯度反演算法(DCG)。把模型更新迭代當(dāng)作是解方程Ax=B,把求解方程轉(zhuǎn)換為求f(x)最小值，表示為：

(11)

把數(shù)據(jù)空間的反演模型迭代形式進(jìn)行分解：

(12)

運(yùn)用共軛梯度算法求解上式方程組，以下為DCG反演計(jì)算流程；

1)選擇較好的mprior和初始模型m0;

2)外層循環(huán)：n=0，1，2…….,計(jì)算dn=d-F(mn)+J(mn-mprior)。F(mn)為當(dāng)前迭代中的模型響應(yīng);

4)設(shè)x0=0,r0=p0=dn;

5)內(nèi)循環(huán)(共軛梯度迭代);

8)xk+1=xk+αkpk;

9)rk+1=rk-αkyk;

10)若|(rk+1)|

13)dn、bn、xk、yk、pk、rk等是數(shù)據(jù)空間共軛梯度反演中的向量，αk和βk是標(biāo)量步長(zhǎng)。在進(jìn)行DCG反演的時(shí)候，對(duì)于Ck的選取非常重要[2]。其中γk是一個(gè)特定的標(biāo)量。選擇合適的Ck，可使共軛梯度算法的迭代次數(shù)大大減少，從而降低反演時(shí)間、提高效率。

2.3 非線性共軛梯度法 (NLCG)反演

隨著大地電磁反演方法的不斷發(fā)展，相關(guān)學(xué)者提出利用CG求解線性方程組，并將該方法用于求解無(wú)約束最優(yōu)化問(wèn)題。非線性共軛梯度法(NonlinearConjugateGradient，NLCG)與CG算法相比，它們的基本理論大體一致。NLCG廣泛應(yīng)用于大地電磁反演中，是解決非線性最優(yōu)化問(wèn)題的一種方法，并已在實(shí)際大地電磁資料解釋中取得了較好的應(yīng)用效果。

NLCG目標(biāo)函數(shù)為

φ(m)=(d-F[m])TV-1(d-F[m])+λmTLTLm

(13)

其中，F(xiàn)為正演函數(shù);V為誤差向量方差的正定矩陣;L是二階差分拉普拉斯算子。

迭代過(guò)程為

(14)

式中，步長(zhǎng)αk利用一維搜索確定，搜索方向經(jīng)過(guò)迭代產(chǎn)生，表達(dá)式為

(15)

(16)

由此可知，NLCG的主要計(jì)算量在于計(jì)算梯度和步長(zhǎng)，步長(zhǎng)的選取標(biāo)準(zhǔn)是使得目標(biāo)函數(shù)滿足一定條件時(shí)就搜索[5]。當(dāng)沒(méi)有先驗(yàn)條件時(shí)， 令Ck=I。

在NLCG算法中，在進(jìn)行每一次迭代時(shí)有下式

(17)

(18)

本文采用另外一種非標(biāo)準(zhǔn)的計(jì)算策略，具體表達(dá)式為：

(19)

該表達(dá)式滿足Wolfe條件，也能夠快速找到局部最小值點(diǎn)，這是它最大的好處。不難得出，在每次迭代中需要計(jì)算梯度一次以及目標(biāo)函數(shù)若干次，這與最速下降法是非常相似的。

表1 NLCG算法Table 1 The flow of NLCG algorithm

2.3 有限內(nèi)存擬牛頓法 (LBFGS)反演

隨著大地電磁反演的發(fā)展，相關(guān)學(xué)者用DFP 算法來(lái)構(gòu)建該近似矩陣，提出一種新型的擬牛頓算法，簡(jiǎn)稱BFGS法，但該方法仍然需存儲(chǔ)矩陣，使得反演時(shí)間增加，效率降低。針對(duì)這一問(wèn)題，相關(guān)學(xué)者在此基礎(chǔ)上，提出L-BFGS法，該算法的優(yōu)點(diǎn)是：構(gòu)建近似的Hessian矩陣，只利用最近一次迭代的曲率信息。

(20)

(21)

通過(guò)重復(fù)計(jì)算式(2)～式(20)式，L-BFGS 近似 Hessian 矩陣的表達(dá)式可表示為[5]：

(26)

(27)

表2 L-BFGS雙遞歸算法流程Table 2 The flow of L-BFGS double recursive algorithm

表3 L-BFGS算法流程Table 3 The flow of L-BFGS algorithm

3 理論模型反演對(duì)比分析

為了驗(yàn)證DCG、NLCG、L-BFGS這三種不同大地電磁二維反演方法各自的優(yōu)缺點(diǎn)，本節(jié)以理論模型試驗(yàn)為例進(jìn)行分析。通過(guò)設(shè)計(jì)兩個(gè)理論地電模型——高、低阻體組合模型和復(fù)雜組合模型，用上述三種方法分別進(jìn)行反演，最后對(duì)比分析反演效果。

3.1 高、低阻體組合模型

設(shè)計(jì)了一個(gè)含有高阻異常體和低阻異常體的組合模型，在模型中圍巖電阻率為100 Ω·m，范圍是均勻半空間。包含兩個(gè)異常體(圖1)，分別是：電阻率為10 Ω·m的低阻異常體，寬度為30 000 m，厚度為15 000 m，水平方向上的位置為-30 000～0 m，深度方向上的位置為0～15 000 m；電阻率為1 000 Ω·m的高阻異常體，寬度為30 000 m，厚度為15 000 m，水平方向上的位置為0～30 000 m，深度方向上的位置為0～15 000 m，如圖1所示。此模型反演所用的數(shù)據(jù)是TE、TM、Tzy三種不同模式下的響應(yīng)數(shù)據(jù)，并且對(duì)電阻率加入了5%的噪音，在相位上加入了1.45度的噪音。其中測(cè)點(diǎn)總共40個(gè)，點(diǎn)距是2 500 m，每個(gè)測(cè)點(diǎn)上的頻點(diǎn)數(shù)為12，頻率范圍為0.001～3.333 Hz。

3.1.1 反演模型結(jié)果對(duì)比分析

反演所用數(shù)據(jù)是，TE、TM、Tzy三種不同模式下的響應(yīng)數(shù)據(jù)，下面一一列舉，見(jiàn)圖2～圖4。

圖2 TE模式下不同方法反演的電阻率斷面Fig.2 Resistivity profile of different inversion methods in TE mode

圖3 TM模式下不同方法反演的電阻率斷面Fig.3 Resistivity profiles retrieved by different methods in TM mode

圖4 Tzy模式下不同方法反演的電阻率斷面Fig.4 Resistivity profiles retrieved by different methods in Tzy mode

對(duì)不同極化模式下的響應(yīng)數(shù)據(jù)采取三種不同方法進(jìn)行反演，得到以上結(jié)果。與真實(shí)模型相比，TM模式的三種反演方法的反演結(jié)果基本上恢復(fù)了異常體的位置、尺寸，低阻異常體的反演電阻率與真實(shí)的非常相近，三種方法反映高、低阻異常體的電阻率值都會(huì)有一定的偏差，異常都會(huì)向四周發(fā)散，不同反演方法都能較好的恢復(fù)地電模型。NLCG和L—BFGS對(duì)高阻異常體的反演較好，反演效果差不多，能大致恢復(fù)高阻異常體的尺寸與電阻率值，但三種方法對(duì)高阻異常體的反演電阻率都偏小。

3.1.2 不同反演算法迭代情況

由表4可知，迭代過(guò)程中，隨著迭代次數(shù)的增加，RMS值都呈下降狀態(tài)，RMS的不斷減小，進(jìn)一步說(shuō)明三種方法收斂性都很好。NLCG比DCG、L-BFGS的RMS值減小更快、收斂速度更快。但是隨著迭代次數(shù)增加，三種反演方法不斷收斂的過(guò)程中，它們的RMS值都趨于1，表示三種反演方法都有較好的收斂性。

表4 高、低阻模型模式反演算法迭代情況

3.2 復(fù)雜模型

設(shè)計(jì)了一個(gè)復(fù)雜模型(圖5)，在模型中研究區(qū)域?yàn)?00 Ω·m的均勻半空間，包含四個(gè)異常體，分別是：

1)電阻率為10 Ω·m的低阻異常體，寬度為30 000 m，厚度為25 000 m，水平方向上的位置為-50 000～-20 000 m，深度方向上的位置為3 000～28 000 m；

2)電阻率為10 Ω·m的低阻異常體，寬度為25 000 m，厚度為15 000 m，水平方向上的位置為-20 000～5 000 m，深度方向上的位置為13 000 m～28 000 m；

3)電阻率為1 000 Ω·m的高阻異常體，寬度為10 000 m，厚度為10 000 m，水平方向上的位置為25 000～35 000 m，深度方向上的位置為4 000～14 000 m；

4)電阻率為1 000 Ω·m的高阻異常體，寬度為10 000 m，厚度為1 000 m，水平方向上的位置為-30 000～-20 000 m，深度方向上的位置為0～1 000 m。

此模型反演所用的數(shù)據(jù)是TE、TM、Tzy三種不同模式下的響應(yīng)數(shù)據(jù)，并且對(duì)電阻率加入了5 %的噪音，在相位上加入了1.45度的噪音。其中測(cè)點(diǎn)總共40個(gè)，點(diǎn)距是2 500 m，每個(gè)測(cè)點(diǎn)上的頻點(diǎn)數(shù)為12，頻率范圍為0.001～3.333 Hz。

3.2.1 反演模型結(jié)果對(duì)比分析

反演所用數(shù)據(jù)為TE、TM、Tzy三種不同模式下的響應(yīng)數(shù)據(jù)，下面一一列舉，見(jiàn)圖6～圖8。

圖6 TE模式下不同方法反演的電阻率斷面Fig.6 Resistivity profile of different inversion methods in TE mode

圖7 TM模式下不同方法反演的電阻率斷面Fig.7 Resistivity profile of different inversion methods in TM mode

圖8 Tzy模式下不同方法反演的電阻率斷面Fig.8 Resistivity profile of different inversion methods in Tzy mode

對(duì)不同極化模式下的響應(yīng)數(shù)據(jù)采取三種不同方法進(jìn)行反演，得到以上結(jié)果。與真實(shí)模型相比，三種反演方法的反演結(jié)果對(duì)異常體的位置、尺寸都可以大致地恢復(fù)出來(lái)，不同反演方法都能較好地恢復(fù)地電模型。Tzy模式時(shí)，相對(duì)于DCG的反演結(jié)果，NLCG和L-BFGS反演結(jié)果更能反映出異常體的尺寸與電阻率的真實(shí)值。在反演的電阻率數(shù)值上，三種方法反映高、低阻異常體的電阻率值都會(huì)有一定的偏差，異常都會(huì)向四周發(fā)散。NLCG和L-BFGS對(duì)高阻異常體的反映較好，能大致確定高阻異常體的尺寸與電阻率值，但三種方法對(duì)高阻異常體的反演電阻率都偏小。L-BFGS和NLCG總體反演效果相當(dāng)，L-BFGS更適合復(fù)雜介質(zhì)模型。

3.2.2 不同反演算法迭代情況

由表5可知，迭代過(guò)程中，隨著迭代次數(shù)的增加，RMS擬合差值都呈下降狀態(tài)，RMS不斷減小，進(jìn)一步說(shuō)明三種方法收斂性都很好。NLCG比DCG、L-BFGS的RMS擬合差值減小更快、收斂速度更快。而隨著迭代次數(shù)增加，三種反演方法不斷收斂的過(guò)程中，它們的RMS值都趨于1，表示三種反演方法都有較好的收斂性。但對(duì)于復(fù)雜地電模型，L-BFGS能獲得更小的RMS擬合差。

表5 復(fù)雜模型反演算法迭代情況

4 結(jié) 論

本文通過(guò)對(duì)大地電磁不同二維反演算法的研究，從麥克斯韋方程組出發(fā)，利用理論模型正演響應(yīng)的數(shù)據(jù)進(jìn)行反演，對(duì)DCG、NLCG、L-BFGS三種反演算法的反演效果進(jìn)行了比較，得到了以下結(jié)論：

1)通過(guò)對(duì)高、低阻模型和復(fù)雜模型的反演對(duì)比可以看出，三種方法均能較好地反演出異常體的位置以及大小，但是對(duì)異常體電阻率的反演值與理論真實(shí)值之間總存在一定偏差，這是不可避免的；同時(shí)，三種方法反演效果大體相當(dāng)，無(wú)論是低阻異常體還是高阻異常體都會(huì)向四周發(fā)散，L-BFGS和NLCG總體反演效果相當(dāng)。

2)通過(guò)對(duì)高、低阻模型和復(fù)雜模型的反演，三種方法在迭代過(guò)程中，數(shù)據(jù)擬合差RMS值都呈不斷下降狀態(tài)。NLCG比DCG、L-BFGS的RMS值減小更快、收斂速度更快。隨著迭代次數(shù)增加，它們的RMS值最后都趨于1，表示三種反演方法都有較好地收斂性，但是對(duì)于復(fù)雜地電模型，L-BFGS的反演效果更好。

3)由于不同反演方法各有自己的特點(diǎn)，所以在實(shí)際選擇方法時(shí)要考慮其優(yōu)缺點(diǎn)，從而更好的發(fā)揮不同反演方法的優(yōu)點(diǎn)；同時(shí)，反演算法的運(yùn)行速度和所需內(nèi)存不一樣，進(jìn)行算法設(shè)計(jì)時(shí)候，應(yīng)該充分考慮到這一點(diǎn)，從而讓大地電磁反演工作效率更高。