周君君， 胡祥云，肖調(diào)杰，白寧波

(1.中國地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，湖北 武漢 430074；2.國防科技大學(xué) 并行與分布處理實(shí)驗(yàn)室，湖南 長沙 410073；3.河南理工大學(xué)物理與電子信息學(xué)院，河南 焦作 454000)

0 引 言

大地電磁測深法(Magnetotelluric，MT)是一種以天然電磁場作為場源的電磁勘探方法，被廣泛地應(yīng)用于固體礦產(chǎn)勘查、石油天然氣勘探等領(lǐng)域[1]。大地電磁勘探資料反演技術(shù)難度大，其處理結(jié)果準(zhǔn)確與否直接影響到地質(zhì)解釋的精度，是MT勘探的關(guān)鍵技術(shù)。大地電磁反演的最終目的是實(shí)現(xiàn)目標(biāo)泛函的最優(yōu)化。反演的方法分為線性反演[2-3 ]和非線性反演[4-6]，兩類算法的區(qū)別在于是否對目標(biāo)函數(shù)進(jìn)行線性化處理。兩類算法各自有自己的缺陷，線性反演算法對初始模型比較依賴，容易陷入局部極小；非線性反演算法對初始模型依賴大幅減弱，但是隨著模型參數(shù)的增加導(dǎo)致搜索空間呈幾何級數(shù)增加。

大地電磁反問題是不適定問題，正則化處理[7]是獲得反問題穩(wěn)定解的有效方法。地球物理工作者將正則化的思想與大地電磁反演問題相結(jié)合，對數(shù)據(jù)目標(biāo)泛函添加不同的穩(wěn)定泛函，如最小模型穩(wěn)定泛函、最光滑模型穩(wěn)定泛函、最小支撐泛函、最小梯度支撐泛函等。

在經(jīng)典的Tikhonov正則化框架下，可以選擇零階、一階、二階差分算子來獲得最小、平緩、光滑解。如基于最平緩模型約束如Occam反演[8-11]、基于最光滑模型約束的NLCG反演[12]，2種反演方法能得穩(wěn)定平滑的反演結(jié)果；然而考慮到對地電界面進(jìn)行清晰的成像，聚焦反演在MT資料處理中也有較多應(yīng)用，能夠獲得較清晰的電性分界面[13-16]。光滑反演與聚焦反演均為正則化的處理方式，兩者的區(qū)別在于正則化項(xiàng)。

大地電磁一維反演正則化反演中，比較典型的有Occam反演[8]和ARIA一維反演[17]。針對大地電磁反演問題的不適定性，特別對于二維和三維的反演，算法的收斂速度占重要地位[18-22]。2種正則化反演算法被提出，最小二乘光滑約束反演和同倫正則反演算法。其中同倫正則算法是將同倫延拓思想與Tikhonov正則化思想相結(jié)合，使算法具備了同倫算法對非線性問題的收斂性好的優(yōu)點(diǎn)。在詳細(xì)敘述2種算法的工作原理之后，采用理論模型和實(shí)測數(shù)據(jù)對2種算法的正確性和可靠性進(jìn)行驗(yàn)證。

1 反演理論和方法

1.1 大地電磁反演理論

反演問題可表示為

d=F(m)

(1)

式中：d為N維觀測數(shù)據(jù)向量；m為M維模型參數(shù)向量；F為正演響應(yīng)算子。

反演問題是已知觀測數(shù)據(jù)向量d，尋求合理的地下模型參數(shù)向量m來擬合觀測數(shù)據(jù)向量d。式(1)是不適定問題，解具有不唯一性和不穩(wěn)定性，因此需要正則化處理獲得反演問題的穩(wěn)定解。根據(jù)Tikhonov理論[7]，構(gòu)造反問題的目標(biāo)泛函如下：

Φ(m)=φ(m)+αs(m)

(2)

式中：Φ(m)為目標(biāo)泛函；φ(m)為數(shù)據(jù)目標(biāo)泛函；s(m)為穩(wěn)定泛函，其作用是使反演過程穩(wěn)定；α為正則因子。

1.2 最小二乘光滑約束反演算法

最小二乘光滑約束反演采用基于先驗(yàn)?zāi)Ｐ偷淖罟饣Ｐ图s束，在L2范數(shù)的框架下，反演的目標(biāo)泛函表示為

(3)

式中：Wd和Wm分別為數(shù)據(jù)和模型加權(quán)矩陣；mapr是先驗(yàn)?zāi)Ｐ汀?/p>

對目標(biāo)泛函式(3)中的F(m)在m-mk處線性化代處理，則式(3)可表示為

Φ(m)≈||Wd(F(mk)+Jk(m-mk)-d)||2+

α||Wm(m-mapr)||2

(4)

式中：Jk為靈敏度矩陣；k為當(dāng)前迭代次數(shù)。

令式(4)的一階變分等于零，得到如下的迭代式,即

1.2 同倫正則反演算法

同倫算法最早用于方程求根，KALLER等[23]首次采用同倫算法求解反問題，提出了一種地震波射線快速追蹤算法。其后，許多學(xué)者對同倫算法進(jìn)行了研究[24-25]。同倫算法由于具有收斂性好的優(yōu)點(diǎn)，因而是求解非線性問題的強(qiáng)有力工具。

同倫算法來源于代數(shù)拓?fù)鋵W(xué)同倫算法是代數(shù)拓?fù)鋵W(xué)中的一個(gè)基本概念，設(shè)均為拓?fù)淇臻g，f和g是X到Y(jié)的連續(xù)映射，若存在連續(xù)映射H:X×[0,1]→Y，使得對任意的x∈X，有H(x,0)=f(x)，H(x,1)=g(x)，則稱f和g同倫，并稱H是連接f(x)和g(x)的一個(gè)同倫。

針對反問題(1)，考慮用不動(dòng)點(diǎn)同倫進(jìn)行求解。將同倫技巧和Tikhonov正則化算法結(jié)合，構(gòu)造了大地電磁反演目標(biāo)泛，該泛函如下：

(5)

式中：m0是初始模型；α正則因子取0<α≤1。

將F(m)在mk處進(jìn)行一階泰勒展開，并令(5)式的一階變分等于零，得到如下的迭代格式，即

mk+1=mk+Δα

(6)

(7)

當(dāng)α=0時(shí)，式(5)為模型目標(biāo)泛函g1(m)=||Wm(m-m0)||2，有唯一解m=m0；當(dāng)α=1時(shí)，g2(m)=||Wd(d-F(m))||2為數(shù)據(jù)目標(biāo)泛函。g1(m)的零點(diǎn)已知，g2(m)的零點(diǎn)是待求的參數(shù)。因此，反問題的解可以表示為從給定的問題到目標(biāo)問題的一個(gè)連續(xù)逼近的過程。應(yīng)用該方法求解反問題式(1)，算法可以全局收斂到反問題的最優(yōu)解，且在一定程度上改善了線性反演對初值的依賴性。

1.3 正則因子的選取

正則因子是數(shù)據(jù)目標(biāo)泛函和模型目標(biāo)泛函的折衷。正則因子較大，則反演結(jié)果過度強(qiáng)調(diào)模型約束；正則因子太小，表示反演結(jié)果過度注重?cái)M合數(shù)據(jù)，從而壓制了模型約束，這將導(dǎo)致反演結(jié)果容易受到噪聲影響產(chǎn)生虛假異常。因此，合適的正則因子選擇至關(guān)重要。在正則因子選取上，已有的工作有正則因子遞降法[26]，即是從初始的正則因子α0出發(fā)，每次迭代根據(jù)一定的步長遞減正則因子。陳小斌2005年提出的MD和CMD方案[17]，其以每一次迭代的數(shù)據(jù)目標(biāo)泛函和模型目標(biāo)泛函的商作為下一次的迭代正則因子的值。及Zhdanov提出的自適應(yīng)算法[27]。

鑒于Morozov偏差原理在正則因子選取上有著廣泛的研究和應(yīng)用，按照Morozov偏差原理選取正則因子，即所求正則因子α的選取與原始數(shù)據(jù)的觀測誤差δ匹配，即滿足以下

||F(m)-d||2=δ2

(8)

其中，δ是觀測數(shù)據(jù)的誤差水平，構(gòu)造泛函如

ψ(α)=||F(m)-d||2-δ2=0

(9)

對(9)式進(jìn)行一階泰勒展開，則

ψ(α)≈l(α)=ψ(αk)+ψ′(αk)(α-αk)

(10)

令l(α)=0，則第k次迭代的Newton迭代格式為

(11)

2 算例分析

2.1 算例一

假設(shè)有3層K型地電模型具體的參數(shù)：視電阻率ρ1=10 Ω·m,ρ2=200 Ω·m,ρ3=20 Ω·m；深度h1=500 m,h2=2 000 m。反演過程中將Bostick反演結(jié)果作為反演的初始模型，最小二乘光滑約束反演結(jié)果如圖1所示。

圖1中的圖標(biāo)中的1表示最小二乘光滑約束反演算法正則因子由遞降法決定[26]；2表示正則因子由Morozov偏差原理選取。由圖1a可得，采用偏差原理選取正則因子的最小二乘光滑約束反演結(jié)果與實(shí)際地電模型最吻合。和Bostick相比能更準(zhǔn)確反映中間高阻體大小和位置。除此之外，還能看出同一種算法，正則因子的選取方法不同，反演效果也不同。本文偏差原理決定正則因子反演結(jié)果(圖1a的紅色曲線)比用遞降法決定正則因子反演結(jié)果(圖1a中藍(lán)綠色曲線)效果好，更吻合真實(shí)模型。

圖1 K型地電模型反演結(jié)果對比

從圖1b和圖1c中可以看出，在數(shù)據(jù)的擬合程度上，采用偏差原理決定正則因子的最小二乘光滑約束反演比Bostick、采用遞降法決定正則因子的最小二乘光滑約束反演的擬合效果要好?？紤]到實(shí)測數(shù)據(jù)通常帶有一定的誤差，為了能夠模擬真實(shí)的觀測數(shù)據(jù)，驗(yàn)證算法的穩(wěn)定性，對理論數(shù)據(jù)添加5%的隨機(jī)高斯噪聲進(jìn)行反演，反演的結(jié)果如圖2所示。由圖2可知，加入5%的隨機(jī)高斯噪聲后反演的結(jié)果基本能反應(yīng)出高阻體的位置和大小，說明了最小二乘光滑約束反演算法能有著良好的魯棒性。圖3給出了無噪聲和 添加5%高斯隨機(jī)噪聲對應(yīng)的迭代收斂曲線圖。最小二乘光滑約束反演算法用偏差原理、遞降法選取正則因子的收斂曲線分別對應(yīng)圖3中的紅色、青綠色曲線?？芍闷钤頉Q定正則因子對應(yīng)的收斂曲線收斂速度更快，迭代次數(shù)少，說明了本文算法一定程度上提高了反演效率。

圖2 K型地電模型的反演結(jié)果對比(添加5%噪聲)

圖3 迭代收斂曲線

2.2 算例二

假設(shè)有4層KH型地電模型具體的參數(shù)：視電阻率ρ1=30、,ρ2=200、ρ3=10、ρ4=100 Ω·m；深度h1=500 m,h2=1 000 m,h3=1 500 m。為了說明觀測數(shù)據(jù)誤差對反演的影響，對原始模型正演的數(shù)據(jù)分別添加了1%、5%、10%、20%的高斯隨機(jī)噪聲，并分別對不同的噪聲水平進(jìn)行同倫正則反演。算法反演的結(jié)果見表1。

表1 四層地電模型的反演結(jié)果

由表1可知，同倫正則反演能很好的反演地電模型參數(shù)，且有著很好的抗噪能力。圖4—圖6分別給出了反演結(jié)果和模型正演結(jié)果在無噪聲、5%噪聲和20%噪聲條件下的視電阻率和相位擬合效果圖。

圖4 KH型地電模型的反演擬合結(jié)果(無噪聲)

圖5 KH型地電模型的反演擬合結(jié)果(5%噪聲)

圖6 KH型地電模型的反演擬合結(jié)果(20%噪聲)

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

為進(jìn)一步驗(yàn)證同倫正則反演算法的有效性，對實(shí)測數(shù)據(jù)進(jìn)行反演。實(shí)測數(shù)據(jù)來自文獻(xiàn)[28]，是吉林樺皮廠地?zé)崽锏囊粋€(gè)實(shí)測點(diǎn)數(shù)據(jù)，用一個(gè)七層的地電模型對該實(shí)測數(shù)據(jù)進(jìn)行反演，將反演結(jié)果和Bostick反演、及比較常用的線性算法Occam的反演的結(jié)果進(jìn)行了對比。圖7給出了3種算法的視電阻率和相位的擬合曲線。針對視電阻率曲線的擬合結(jié)果來看，雖然Occam反演與同倫正則反演較結(jié)果整體優(yōu)于Bostick算法。但是Occam的高頻部分?jǐn)M合較差，而同倫正則反演算法在高低頻的擬合效果都較好；從相位曲線的擬合結(jié)果來看，Bostick在高低頻的擬合效果較好，中頻的擬合效果較差，Occam在高低頻段擬合效果較差，而同倫正則算法效果是最佳的。說明了本文算法用于實(shí)測大地電磁的反演也十分有效。

圖7 實(shí)測數(shù)據(jù)反演結(jié)果對比

3 結(jié) 論

1)給出了2種正則化算法：最小二乘光滑約束反演與同倫正則反演算法。其中，首次將同倫正則反演應(yīng)用到大地電磁反演中。2種算法應(yīng)用到不同地電模型中，均取得較好的反演結(jié)果。體現(xiàn)了算法的有效性和穩(wěn)健性。

2)正則因子的選取上，兩種算法均采用偏差原理決定正則因子。算例一將其與傳統(tǒng)正則因子遞減法的反演結(jié)果進(jìn)行了對比分析，證明了用偏差原理選取正則因子，效率更高，反演效果更好。

3)同倫正則反演算法是考慮到同倫算法收斂性好的優(yōu)點(diǎn)，構(gòu)造了目標(biāo)泛函，并成功應(yīng)用在大地電磁一維反演中。將來，可以嘗試將此算法推廣到大地電磁二維、三維的反演中。