大地電磁非線性共軛梯度一維反演

管貽亮,李予國,胡祖志,黃江波,胡祥云

(1.中國地質大學(武漢)地球物理與空間信息學院,湖北武漢430074;2.中國海洋大學海洋地球科學學院,山東青島266100;3.中國石油集團東方地球物理勘探有限責任公司,河北涿州072751)

大地電磁測深(MT)是以天然交變電磁場為場源,研究地殼和上地幔構造的一種地球物理勘探方法[1]。20世紀50年代初，蘇聯學者Tikhonov[2]和法國學者Cagniard[3]建立了探測地球深部電性結構的大地電磁測深法。該方法工作簡便,勘探深度大,在礦產普查與勘探、地殼與上地幔電性結構研究、地熱田的調查以及石油天然氣勘探等方面發(fā)揮了重要作用[4]。

當前大地電磁測深一維反演常用的反演方法有Bostick近似反演法、高斯-牛頓法、馬夸特法[5]、模擬退火法[6]、Occam法[7]等,不同的方法有各自不同的特點[8-9]。在高斯-牛頓算法和Mackie-Madden算法[10-11]的基礎上發(fā)展起來的非線性共軛梯度法(NLCG)多用于二維反演[12],Newman等[13]將其用于三維反演。它突破了線性化迭代反演的框架,能夠直接對非二次極小化問題求解[14]。但是該算法由于計算量較大,正則化因子選取較困難,很少被用于實際資料的一維反演,在一維反演方面的應用研究也很少。

本文結合一維反演實際,提出了非線性共軛梯度一維反演算法的簡化和改進方案。

1 大地電磁測深的基本理論

大地電磁測深的場源是天然交變電磁場,其理論基礎是卡尼爾-吉洪諾夫理論,將場源近似看成平面波垂直入射大地。當交變電磁場在地下介質中傳播時,由于趨膚深度的作用,不同頻率的信號具有不同的穿透深度,在地面上觀測大地電磁場,它的頻率響應能反映地下介質電性的垂向分布情況[15]。為了得到地下介質電阻率的分布信息,引入波阻抗的概念:Z=E/H。

地面波阻抗可依次用其下相鄰介質頂面的波阻抗來表示,由于波阻抗在分界面上是連續(xù)的,對于層狀介質,任一層頂界面的波阻抗等于其上相鄰介質底界面的波阻抗,因而問題轉變?yōu)榍蠼馔粚咏橘|頂、底界面波阻抗之間的關系[1]。因為每層介質都是均勻各向同性的,因此頂、底界面的波阻抗可用積分常數聯系起來,由此得到大地電磁正演的遞推公式

(1)

式中：Z為波阻抗;E為電場強度;ω為圓頻率;k為復波數;h為地層厚度;μ為磁導率;ρa為視電阻率。

2 大地電磁數據NLCG一維反演算法

2.1 非線性共軛梯度法原理

Rodi等[12]將NLCG法反演問題表示成

(2)

式中:d為數據矢量；m為模型矢量；F為正演算子；e是殘差向量。

按照Tikhonov和Arsenin的正則化方法[16]處理目標函數的最小化問題,定義目標函數Ψ為

(3)

正演函數的雅可比矩陣A為

(4)

目標函數Ψ(m)的梯度向量g和Hessian矩陣H定義為

(5)

式中：B是F的Hessian矩陣；M是模型參數個數；N是地層個數。在非線性共軛梯度算法中,同樣首先給定初始模型,通過求解一元極小化問題或者沿著計算的方向進行線性搜索確定模型序列：

(6)

式中：m0為給定的初始模型；pl為搜索方向；al為搜索步長。在每次迭代極小化目標泛函時,需要沿搜索方向pl計算搜索步長al,搜索方向的迭代格式與線性共軛梯度近似:

(7)

為了簡化計算,我們選取Polak-Ribiere共軛梯度法[19]來極小化目標函數,則(7)式中有

(8)

其中，Cl為預調節(jié)因子。

在線性搜索過程中,如果搜索收斂,則新的搜索方向pl由(7)式計算;如果搜索失敗,則pl取為最速下降方向(pl=-Clgl),打破和前面搜索方向的共軛性。

2.2 NLCG算法簡化

在二維算法中,Rodi等[12]和Newman等[20]分別給出了預調節(jié)因子的構造方案,其中Rodi等用如下方法構造預調節(jié)因子：

(9)

式中：γl是特定的標量。應用線性共軛梯度法,通過對方程組(10)中h的求解將預調節(jié)因子作用于梯度向量：

(10)

在具體實現時,這需要增加很大的計算量,考慮到一維反演的復雜程度遠低于二維,為了使反演算法更加快速,嘗試在一維反演中將Cl取為單位矩陣,即不再進行預調節(jié)。

最優(yōu)步長是反演中最重要的變量。二維反演中使用單變量高斯-牛頓線性搜索方法,在每一次共軛梯度循環(huán)中進行多次高斯-牛頓迭代,通過與模型進行實時對比確定最優(yōu)步長。這種方法雖然能夠得到較好的步長,但是需要重復調用正演函數及其它反演參數,對一維反演而言增加了太多的計算量,得不償失。因此我們嘗試用單步線性搜索法求步長,即

(11)

這樣求得的步長雖然不是最優(yōu)結果,但經驗證其反演效果能夠達到預期,最主要的是減少了很多的計算量,相對于高斯-牛頓迭代節(jié)省了至少1/3的時間。就一維反演而言這樣的簡化是合適的,但對于二維以上大數據量的反演,這樣的簡化會減緩目標函數的收斂速度,增加迭代次數,影響反演效率,甚至導致反演失敗。

簡化后的NLCG算法在每次迭代過程中采用Polak-Ribiere共軛梯度法計算線性搜索方向,在每一搜索方向上用單步線性搜索法求取最優(yōu)步長,得到新的模型,在規(guī)定擬合差或迭代次數范圍內結束循環(huán),詳細流程如圖1所示。

圖1 非線性共軛梯度反演流程

2.3 理論模型算例分析

設計了多個理論模型來驗證NLCG反演算法的有效性。以7層水平層狀介質模型為例,將正演視電阻率計算值加5%的隨機噪聲作為觀測值,初始模型層厚按一定比例遞增,每層電阻率取視電阻率平均值,根據反演結果不斷調整正則化因子,以使反演效果達到最佳。同時,為避免模型出現較大的跳躍,嘗試對增量進行壓制(增量乘以限制系數),這樣得到的反演曲線比較光滑。

圖2a為7層模型反演結果，可以看出,NLCG算法的適應性和穩(wěn)定性比較好,在有一定噪聲干擾的情況下,電性層的劃分以及每層的電阻率值均能較為準確地與設計模型對應,尤其是淺層很薄的電性層依舊能被準確識別,且不會出現冗余構造,說明NLCG算法反演精度和分辨率較高,在小構造的識別中具有一定的優(yōu)勢。圖2b中擬合曲線能夠較好反映觀測數據中微小的變化,說明正則化因子的存在很好地兼顧了模型粗糙度和擬合度。

圖2 7層水平層狀介質模型反演結果a 反演結果; b 擬合曲線

圖3為3層和4層層狀介質模型的目標函數收斂曲線,可以看出,對于不同的模型,大約迭代10次以后目標函數就能夠趨于收斂且穩(wěn)定,說明NLCG算法的收斂速度較快,所需的迭代次數較少,在計算方面具有一定的優(yōu)勢。

NLCG算法能夠節(jié)省運算時間、提高反演效率,同時又能保證反演的分辨率,因此適用于大規(guī)模的非線性一維反演計算。

圖3 目標函數收斂曲線

3 NLCG一維反演算法改進

通過模型算例分析可以看出,NLCG算法反演效果很好,具有較高的分辨率。但是該算法需要滿足一些前提條件,即最優(yōu)的初始模型和最合適的正則化因子,尤其是簡化后的算法對初始模型的依賴度很高。為了達到最好的反演效果,在每次試驗中都需要不斷調整初始模型和正則化因子,這在一定程度上降低了反演的效率,增加了工作量。因此,以最小的計算量和最短的時間來選擇最優(yōu)的初始模型和正則化因子成為提高NLCG算法反演效率的關鍵。

3.1 基于Bostick的NLCG反演

由于NLCG算法對初始模型的依賴度很高,因此需要找到一種簡便的反演算法為其提供初始模型。Bostick反演法無疑是最好的選擇。

Bostick反演作為一種近似的反演方法,使用相位進行反演，基本公式為

(12)

式中：H為地層深度,φ為相位。根據測得的視電阻率和相位,用(12)式可以計算出深度-電阻率模型。該方法計算速度快,不需要初始模型,也不需要反復迭代[21],雖然反演精度很低[22],卻是一種為其它反演方法提供初始模型的快捷方法。

3.2 正則化因子的自適應迭代調整

正則化因子作為權重因子起到了平衡目標函數和模型約束的作用。當它比較大時,模型比較光滑,但數據擬合度較差,較小時主要擬合實測數據,模型比較粗糙。因此選擇合適的正則化因子成為提高反演效果的關鍵。

在上述算法中,正則化因子是在程序中給定的,按照反演效果不斷進行修改,這樣需要大量的運算,降低了算法的效率。因此,我們嘗試自動調整正則化因子，以加快算法的運行效率。

考慮到正則化因子在極小化數據擬合差范數和模型范數之間起到折中的作用,參考陳小斌等[23]提出的自適應正則化因子調整方法,結合算法實際,提出了如下調整方案。

令公式(3)中Ψ1=eTV-1e;Ψ2=mTLTLm,則Ψ=Ψ1+λΨ2。其中Ψ1為觀測數據目標函數;Ψ2為先驗約束條件的目標函數,即模型粗糙度函數;λ為正則化因子。以迭代的方式調整正則化因子：

(13)

從公式(13)可以看出,雖然在每次共軛梯度循環(huán)時重新計算了正則化因子,但是不會增加新的計算量,因為計算所需的Ψ1和Ψ2在共軛循環(huán)中已經計算過。

設計多個理論模型驗證了改進算法的有效性,同時驗證了需要給定的正則化因子初始值λ0。以4層KH模型為例,經過不斷試驗,發(fā)現大多數的模型只需將初始值給定1～5中的一個數即可。

圖4為4層KH模型改進后NLCG算法反演結果?？梢钥闯?改進后的NLCG算法具有較好的分辨率,反演結果能夠與模型基本對應,分層效果較好,尤其是低阻層精度比較高。

改進后的NLCG算法不再需要人為地對初始模型和正則化因子進行調整,加快了反演的速度,但是該算法對于高阻層反演精度相對較低,而且由于Bostick反演模型對于深層介質厚度變化反應較快,導致深層電性層劃分具有一定的偏差。

3.3 反演優(yōu)化

NLCG反演對于高阻層的分辨率較差的情況具體表現為淺層高阻出現震蕩,深層高阻擬合精度差(圖2,圖4)。增加迭代次數可以使深層高阻反演精度有所提高,但是淺層高阻震蕩加劇。迭代次數過少則深層高阻反演精度較差,二者不可兼顧,且隨著模型復雜度增加,高阻層的分辨率會進一步降低。所以，我們對反演結果做一次優(yōu)化,對優(yōu)化后的結果再做一次反演(圖5)。

圖4 4層KH模型改進后NLCG算法反演a Bostick反演的初始模型; b 改進后的NLCG反演結果; c 擬合曲線

圖5 NLCG反演優(yōu)化效果a NLCG反演結果; b 二階微分曲線; c 優(yōu)化NLCG反演結果

先做5～10次的NLCG反演,得到較為粗糙的反演模型,對反演曲線做擬合后求微分(二次微分),利用微分曲線極值點(零點)對應反演曲線拐點劃分出地下介質的電性層，確定電性分界面；然后將每一電性層繼續(xù)分成若干薄層,低阻層的劃分可以較厚,高阻層的劃分相對較薄,薄的高阻層劃分有利于快速收斂,高、低阻層差異性劃分也能保證電性分界面的清晰，由于電性層的劃分不完全準確,相鄰電性層的厚度二次劃分不宜差距過大，每一層的電阻率可以取該層反演電阻率平均值,也可統(tǒng)一取視電阻率平均值。將優(yōu)化后的結果作為初始模型再進行一次反演,只需較少的迭代次數就可以達到很高的反演精度,使反演結果的分辨率進一步提高(圖5)。

4 實際資料處理效果分析

為了驗證NLCG算法對于實測資料反演的有效性,我們對A地區(qū)井旁測點的大地電磁數據進行了反演,并與鉆井資料和電測井資料進行了對比。以Bostick反演結果為初始模型,未進行二次反演優(yōu)化，結果如圖6所示。

圖6 A地區(qū)實際MT資料反演效果分析a 實際資料反演結果; b 鉆井巖性及電測井曲線

從反演結果與鉆井數據對比可以看出,反演結果與地下地層結構基本對應：第1層為高阻層,對應于砂巖;第2層為低阻層,反映的是泥巖;第3層對應高阻花崗巖,電阻率比較高。實際資料反演結果表明,改進的NLCG反演算法是有效且比較精確的,反演速度較快,與理論模型反演得出的結論一致。將改進后的NLCG算法用于實際資料電性結構的劃分以及擬二維反演解釋，結果是可靠的。

5 與其它算法的比較

傳統(tǒng)的反演方法中,梯度法(最速下降法)只用到模型的梯度,即一階導數,算法過于簡單,在極小點附近收斂不好;高斯-牛頓算法雖然在極小點收斂比較快,但是其計算量與模型參數量和觀測數據量成正比,在大數據量反演中需要的內存比較大。傳統(tǒng)方法最大的不足是只考慮了擬合優(yōu)度而沒有考慮模型的粗糙度。NLCG算法既保證了在極小點附近收斂比較快,又引入了模型粗糙度,反演的穩(wěn)定性較好,分辨率較高。

以5層層狀介質模型為例,將改進的NLCG算法與馬夸特法[5]進行了對比(圖7)。采用相同的初始模型：層厚按比例遞增,每層電阻率取為視電阻率平均值。由圖7可見,馬夸特法對于深部地層的分辨率較低,NLCG法優(yōu)于馬夸特法,尤其是低阻層精度較高,而且隨著模型復雜度增加,NLCG算法的優(yōu)勢越來越明顯。

采用改進的NLCG法和Occam反演法對A地區(qū)井旁測點的大地電磁數據進行了反演對比,兩種方法都是以Bostick反演結果為初始模型,迭代35次,結果如圖8所示。對比NLCG和Occam反演結果可以看出,二者與觀測數據的擬合都比較好,NLCG反演的結果分層更加清晰,且不會出現Occam反演結果所顯示出的異常高阻層,反演結果更接近實際,在地下電性結構的劃分中可靠度和可信度更高。

圖7 5層水平層狀模型的馬夸特法(a)和NLCG法(b)反演結果對比

圖8 A地區(qū)實際MT資料反演結果對比a Occam與NLCG反演結果; b 擬合曲線

6 結束語

理論研究和實際資料處理結果表明，經過優(yōu)化的NLCG算法能夠節(jié)省運算時間、提高反演效率，同時又能保證反演的分辨率，適合大規(guī)模的非線性一維反演計算。

1) 根據目標函數自動調整正則化因子，避免了對正則化因子的不斷嘗試反演，提高了運行效率。

2) 簡化了預處理因子和最優(yōu)步長的計算，用單步線性搜索法代替迭代法求最優(yōu)步長，使NLCG算法更加高效快捷。

3) 對一次反演結果優(yōu)化后進行二次反演，進一步提高了反演結果的分辨率。

致謝:在與中國海洋大學羅鳴探討算法的過程中,筆者受益匪淺,在此向他表示感謝。

參 考 文 獻

[1] 石英俊,劉國棟,吳廣耀,等.大地電磁測深法教程[M].第1版.北京:地質出版社,1985:1-110

Shi Y J,Liu G D,Wu G Y,et al.The magnetotelluric sounding tutorial[M].1st ed.Beijing:Geological Publishing House,1985:1-110

[2] Tikhonov A N.On determining electrical characteristics of the deep layers of the earth’s crust[J].Dokl Akad Nauk,1950,73(2):95-297

[3] Cagniard L.Basic theory of the magnetotelluric method of geophysics prospecting[J].Geophysics,1953,18(3):605-635

[4] 陳樂壽,王光鍔.大地電磁測深法[M].第1版.北京:地質出版社,1990:1-56

Chen L S,Wang G E.The magnetotelluric sounding method[M].1st ed.Beijing:Geological Publishing House,1990:1-56

[5] Marquardt D W.An algorithm for least-squares estimation of nonlinear parameters[J].Journal of the Society of Industrial and Applied Mathematics,1963,11(2):431-441

[6] 師學明,王家映.一維層狀介質大地電磁模擬退火反演法[J].地球科學,1998,23(5):542-545

Shi X M,Wang J Y.One dimensional magnetotelluric sounding inversion using simulated annealing[J].Earth Science,1998,23(5):542-545

[7] Constable S C,Parker R L,Constable G C.Occam’s inversion:a practical algorithm for generating smooth models from electromagnetic sounding data[J].Geophysics,1987,52(3):289-300

[8] 韓波,胡祥云,何展翔.大地電磁反演方法的數學分類[J].石油地球物理勘探,2012,47(1):177-187

Han B,Hu X Y,He Z X.Mathematical classification of magnetotelluric inversion method[J].Oil Geophysical Prospecting,2012,47(1):177-1877

[9] 陳向斌,呂慶田,張昆.大地電磁測深反演方法現狀與評述[J].地球物理學進展,2011,26(5):1607-1619

Chen X B,Lv Q T,Zhang K.Review of magnetotelluric data inversion methods[J].Progress in Geophysics,2011,26(5):1607-1619

[10] Madden T M,Mackie R L.Three-dimensional magnetotelluric modeling and inversion[J].Proceedings of the IEEE,1989,77(2):318-333

[11] Mackie R L,Madden T R.Three-dimensional magnetotelluric inversion using conjugate gradients[J].Geophysical Journal International,1993,115(1):215-229

[12] Rodi W,Mackie R L.Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion[J].Geophysics,2001,66(1):174-187

[13] Newman G A,Alumbaugh D L.Three-dimensional magnetotelluric inversion using non-linear conjugate gradients[J].Geophysical Journal International,2000,140(2):410-424

[14] 李衛(wèi)東.大地電磁反演方法評析及其在地熱勘探中的應用研究[D].北京:中國地質大學,2006

Li W D.Assessment of magnetotelluric sounding inversion method and research for geothermal exploration[D].Beijing:China University of Geosciences,2006

[15] 柳建新,童孝忠,郭榮文,等.大地電磁測深法勘探—資料處理、反演與解釋[M].第1版.北京:科學出版社,2012:1-111

Liu J X,Tong X Z,Guo R W,et al.The magnetotelluric sounding exploration—data processing,inversion and interpretation[M].1st ed.Beijing:Science Press,2012:1-111

[16] 吉洪諾夫.不適定問題的解法[M].王秉忱譯.第1版.北京:地質出版社,1979:1-160

Tikhonov A N.Solution methods for the ill-posed problems[M].Wang B C translator.1st ed.Beijing:Geological Publishing House,1979:1-160

[17] 胡祖志,胡祥云,何展翔.大地電磁非線性共軛梯度擬三維反演[J].地球物理學報,2006,49(4):1226-1234

Hu Z Z,Hu X Y,He Z X.Pseudo-three-dimensional magnetotelluric inversion using nonlinear conjugate gradients[J].Chinese Journal of Geophysics,2006,49(4):1226-1234

[18] 萬漢平.大地電磁測深TE和TM極化模式對比研究[D].四川:成都理工大學,2010

Wan H P.Comparison of magnetotelluric TE and TM polarization modes[D].Sichuan:Chengdu University of Technology,2010

[19] 王家映.地球物理反演理論[M].第2版.北京:高等教育出版社,2002:1-182

Wang J Y.Inverse theory in geophysics[M].2nd ed.Beijing:Higher Education Press,2002:1-182

[20] Newman G A,Boggs P T.Solution accelerators for large-scale three-dimensional electromagnetic inverse problems[J].Inverse Problems,2004,20(6):151-170

[21] 周虬.一種簡易的一維大地電磁測深反演方法—博斯蒂克法反演及其應用[J].石油地球物理勘探,1985,20(1):80-88

Zhou Q.A simple one-dimensional magnetotelluric inversion methods—Bostick inversion and its application[J].Oil Geophysical Prospecting,1985,20(1):80-88

[22] 徐新學,趙憲堂,王寶勛,等.大地電磁測深資料常規(guī)反演方法的應用現狀及進展[J].地質調查與研究,2004,27(增刊):71-74

Xu X X.Zhao X T,Wang B X,et al.The application status and progress of magnetotelluric sounding data regular inversion[J].Geological Survey and Research,2004,27(S):71-74

[23] 陳小斌,趙國澤,湯吉,等.大地電磁自適應正則化反演算法[J].地球物理學報,2005,48(4):937-945

Chen X B,Zhao G Z,Tang J,et al.An adaptive regularized inversion algorithm for magnetotelluric data[J].Chinese Journal of Geophysics,2005,48(4):937-945