大地電磁一維正則化反演算法研究

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

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

0 引 言

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

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

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

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

1 反演理論和方法

1.1 大地電磁反演理論

反演問題可表示為

d=F(m)

(1)

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

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

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

(2)

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

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

最小二乘光滑約束反演采用基于先驗模型的最光滑模型約束，在L2范數(shù)的框架下，反演的目標泛函表示為

(3)

式中：Wd和Wm分別為數(shù)據(jù)和模型加權矩陣；mapr是先驗模型。

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

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

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

(4)

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

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

1.2 同倫正則反演算法

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

同倫算法來源于代數(shù)拓撲學同倫算法是代數(shù)拓撲學中的一個基本概念，設均為拓撲空間，f和g是X到Y的連續(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)的一個同倫。

針對反問題(1)，考慮用不動點同倫進行求解。將同倫技巧和Tikhonov正則化算法結合，構造了大地電磁反演目標泛，該泛函如下：

(5)

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

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

mk+1=mk+Δα

(6)

(7)

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

1.3 正則因子的選取

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

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

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

(8)

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

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

(9)

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

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

(10)

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

(11)

2 算例分析

2.1 算例一

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

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

圖1 K型地電模型反演結果對比

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

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

圖3 迭代收斂曲線

2.2 算例二

假設有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%的高斯隨機噪聲，并分別對不同的噪聲水平進行同倫正則反演。算法反演的結果見表1。

表1 四層地電模型的反演結果

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

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

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

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

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

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

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

3 結 論

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

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

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