一一種改進(jìn)的二維MT預(yù)條件非線性共軛梯度反演方法

相 鵬

(勝利油田西部新區(qū)研究院,山東東營(yíng) 257000)

一一種改進(jìn)的二維MT預(yù)條件非線性共軛梯度反演方法

相 鵬

(勝利油田西部新區(qū)研究院,山東東營(yíng) 257000)

在大地電磁反演方法中反演精度與計(jì)算效率問題是一對(duì)矛盾,高斯牛頓類方法反演精度高但計(jì)算效率低,非線性共軛梯度類方法計(jì)算效率高,但是反演精度不如高斯牛頓法高。在前人研究的基礎(chǔ)上,提出一種改進(jìn)的預(yù)條件非線性共軛梯度法,通過構(gòu)建性狀更接近高斯牛頓Hessian矩陣的預(yù)條件算子提高反演精度和計(jì)算速度。同時(shí)采用正則化參數(shù)的自適應(yīng)更新算法保證反演穩(wěn)定性和反演精度的平衡。模型實(shí)驗(yàn)驗(yàn)證了該方法的正確性。與其他方法的對(duì)比結(jié)果表明,該方法在保證反演精度的同時(shí),提高了計(jì)算效率。對(duì)中國(guó)西部某地的實(shí)測(cè)MT數(shù)據(jù)進(jìn)行處理解釋的結(jié)果表明,該方法在解決復(fù)雜構(gòu)造問題方面具有較高的實(shí)用價(jià)值。

大地電磁;預(yù)條件非線性共軛梯度;正則化參數(shù);自適應(yīng)算法;反演

近年來,大地電磁(MT)反演技術(shù)隨著數(shù)值計(jì)算技術(shù)的發(fā)展已經(jīng)從一維地電結(jié)構(gòu)反演發(fā)展到二維、三維反演階段。Rodi和Mackie[1]開發(fā)了非線性共軛梯度法(nonlinear conjugate gradient,簡(jiǎn)稱NLCG), NLCG用單位對(duì)角陣作為預(yù)條件算子,但該預(yù)條件算子完全損失了Hessian矩陣信息。Ha和Shin[2]通過截?cái)嘣糎essian矩陣,并在每行保留大量非零元素來構(gòu)建近似Hessian矩陣。這種近似Hessian矩陣用于計(jì)算高斯牛頓法的搜索方向,但是結(jié)果被假象污染,而且計(jì)算效率改善程度有限。筆者提出一種基于稀疏預(yù)條件算子的非線性共軛梯度方法(sparse preconditioned NLCG,簡(jiǎn)稱SP-NLCG)。SPNLCG在近似Hessian矩陣中使用非常少的非零元素,可以保持與NLCG法幾乎相同的計(jì)算開銷。另外,前人研究中鮮見對(duì)正則化參數(shù)的討論,筆者提出一種正則化參數(shù)的自適應(yīng)更新算法,以保證不適應(yīng)條件下反演的穩(wěn)定性及最大限度提高反演結(jié)果的分辨率。

1 正演計(jì)算

對(duì)于二維地電模型,取走向?yàn)閤軸,y軸與x軸垂直,z軸垂直向下(圖1),介質(zhì)模型的電性參數(shù)隨y軸和z軸都發(fā)生變化,而沿走向x軸的電性參數(shù)不發(fā)生變化,即δ δx=0。

由于平面電磁波以任何角度入射地面都以平面波形式垂直地向下傳播,將麥克斯韋方程組以x分量為準(zhǔn)可以得到TM模式和TE模式兩個(gè)獨(dú)立的方程組。

式中,E為電場(chǎng)強(qiáng)度;H為磁場(chǎng)強(qiáng)度;σ為電導(dǎo)率;μ為磁導(dǎo)率;ε為介電常數(shù)。在二維大地電磁正演計(jì)算時(shí),邊界條件是關(guān)系到解的精度的關(guān)鍵,TE模式外邊界條件定義如下:

由于上邊界AB離地面取得足夠遠(yuǎn),使異常場(chǎng)u在AB上為零,所以該處的u=1:

由于左右邊界AD、BC與局部不均勻體距離取得足夠遠(yuǎn),電磁場(chǎng)在AD、BC上左右對(duì)稱,其上的邊界條件為

由于下邊界CD以下為均質(zhì)巖石,局部不均勻體的異常場(chǎng)在CD上為零,CD處的邊界條件為

TM模式外邊界條件定義如下:

(1)由于上邊界AB直接取在地面上,該處的u為1:

(2)左右邊界AD、BC的邊界條件等于TE模式的左右邊界條件。

(3)下邊界CD以下的邊界條件等于TE模式的邊界條件。

對(duì)于TE極化模式,復(fù)數(shù)視電阻率公式為

式中,〈Ex〉為觀測(cè)點(diǎn)電場(chǎng)的平均值;〈Hy〉為觀測(cè)點(diǎn)磁場(chǎng)的平均值。

對(duì)于TM極化模式,復(fù)數(shù)視電阻率公式為

麥克斯韋方程可以用有限差分方程近似[3],用有限差分法求解大地電磁正演問題在一定情況下會(huì)產(chǎn)生偽解[4],當(dāng)頻率或電導(dǎo)率趨近于零時(shí),傳導(dǎo)電流場(chǎng)合磁場(chǎng)滿足零散度條件,而偽解的傳導(dǎo)電流場(chǎng)和磁場(chǎng)的散度卻不為零。為了消除偽解,必須在介質(zhì)和空氣分別強(qiáng)制加入散度條件,目前常用的處理辦法是強(qiáng)加散度校正條件:▽·E=0。

2 稀疏預(yù)條件算子非線性共軛梯度法

2.1 目標(biāo)函數(shù)的構(gòu)建

根據(jù)Tikhonov提出的正則化思想[5-6],大地電磁反演問題的目標(biāo)函數(shù)可以用以下形式構(gòu)建:

式中,d為數(shù)據(jù)矢量;m為模型矢量;正實(shí)數(shù)λ為正則化參數(shù);正定矩陣V是數(shù)據(jù)的協(xié)方差矩陣。式(9)右邊第二項(xiàng)是穩(wěn)定因子項(xiàng),很多學(xué)者對(duì)穩(wěn)定因子進(jìn)行了研究,如最光滑泛函[7-8]和最小梯度支撐泛函[9-11],本文中選擇L為二階差分算子實(shí)現(xiàn)最光滑模型約束反演。

2.2 非線性共軛梯度法

目標(biāo)函數(shù)優(yōu)化的常用方法有最小二乘法、高斯牛頓法、最速下降法等,這些方法都屬于“牛頓型”方法,即模型的搜索方向p由下式求得:

式中,H為目標(biāo)函數(shù)的二階導(dǎo)數(shù),為Hessian矩陣;g為目標(biāo)函數(shù)的梯度。牛頓型方法需要計(jì)算Hessian矩陣及其逆矩陣,在求解大規(guī)模反問題時(shí),內(nèi)存需求和計(jì)算成本都非常高。Fletcher和Reeves[12]提出了NLCG法無需計(jì)算Hessian矩陣,只需計(jì)算目標(biāo)函數(shù)梯度,極大地降低了對(duì)計(jì)算機(jī)內(nèi)存的要求。Rodi等將其應(yīng)用到MT二維反演中,并比較了高斯牛頓法、MM算法和NLCG法的性能,實(shí)驗(yàn)結(jié)果表明NLCG法較其他兩種方法計(jì)算速度更快[1]。

NLCG的基本原理是,首先通過計(jì)算該目標(biāo)函數(shù)在某參考模型處的負(fù)梯度來得到共軛梯度,并在這個(gè)方向上而不是簡(jiǎn)單的梯度方向上搜索一個(gè)最優(yōu)步長(zhǎng)使目標(biāo)函數(shù)達(dá)到極小值。這樣得到的解就應(yīng)該是目標(biāo)函數(shù)在目前共軛梯度方向的極點(diǎn)。其迭代計(jì)算過程可參考文獻(xiàn)[1],NLCG與高斯牛頓法相比,最大的區(qū)別在于求取搜索方向時(shí),不需要求取高斯牛頓Hessian矩陣并對(duì)其求逆矩陣,而是用預(yù)條件算子C代替Hessian矩陣,搜索方向的公式為

因此預(yù)條件算子C的選擇對(duì)于反演的精度和收斂速度有著至關(guān)重要的作用,Rodi等給出預(yù)條件算子為Ci=(γiI+λLTL)-1,γi為特定的正實(shí)數(shù),這種預(yù)條件算子的優(yōu)點(diǎn)是實(shí)現(xiàn)簡(jiǎn)單,無需顯式求解逆矩陣,但是它完全忽略掉了Hessian矩陣的信息,需要很多次迭代才能收斂至預(yù)定條件。

2.3 稀疏預(yù)條件算子

SP-NLCG的算法流程與NLCG相同,區(qū)別在于預(yù)條件算子的構(gòu)建策略。預(yù)條件算子C的性狀與高斯牛頓法Hessian矩陣越接近,反演的精度越高,需要的迭代次數(shù)越少。如果令式(12)中的預(yù)條件算子C等于高斯牛頓法Hessian矩陣,且去掉右邊第二項(xiàng),那么步長(zhǎng)變成單位步長(zhǎng),SP-NLCG退化為高斯牛頓法。顯然,如果C為單位矩陣,SP-NLCG變?yōu)镹LCG。因此,SP-NLCG的收斂速率應(yīng)該介于NLCG和高斯牛頓法之間。

構(gòu)建預(yù)條件算子的主要挑戰(zhàn)在于如何設(shè)計(jì)近似Hessian矩陣,使之有效且計(jì)算廉價(jià)。出于效率的考慮,近似Hessian矩陣必須高度稀疏同時(shí)保持精度。為了得到高度稀疏的近似Hessian矩陣,選擇高斯牛頓Hessian矩陣中最主要的元素,舍棄貢獻(xiàn)小的元素。Hessian矩陣的每個(gè)元素是雅可比矩陣兩列的內(nèi)積,而雅可比矩陣的每一列代表觀測(cè)矢量對(duì)離散化的反演區(qū)域中某個(gè)網(wǎng)格的模型參數(shù)的梯度。實(shí)驗(yàn)表明,當(dāng)模型網(wǎng)格在空間上接近,則對(duì)應(yīng)的雅可比矩陣列向量有很強(qiáng)的相關(guān)性。最強(qiáng)的相關(guān)性為自相關(guān)。這些空間接近的網(wǎng)格產(chǎn)生了Hessian矩陣中的主控元素。因此,在構(gòu)建高度稀疏Hessian矩陣的第i行時(shí),只要選擇包圍著對(duì)應(yīng)雅可比矩陣第i列的網(wǎng)格的那些網(wǎng)格的雅可比列向量即可。選擇的雅可比列向量不必是鄰近的,雖然它們對(duì)應(yīng)的網(wǎng)格在空間上是鄰近的。

稀疏近似Hessian矩陣的結(jié)構(gòu)如圖2(a)所示。

對(duì)角線元素是雅可比矩陣列向量的自相關(guān),對(duì)應(yīng)圓圈區(qū)域的中心網(wǎng)格,而副對(duì)角線元素是圓形區(qū)域中不同網(wǎng)格的互相關(guān)。稱圖2(b)的圓形為Hessian圓,該圓僅包括4個(gè)相鄰網(wǎng)格(二維), Hessian圓越大,近似Hessian矩陣的稀疏性越差。數(shù)值模擬表明在近似Hessian矩陣中包含更多非零項(xiàng)超過一定界限時(shí),不會(huì)使高斯牛頓Hessian矩陣的近似逆矩陣的精度提高更多。實(shí)際上,在多數(shù)情況下精度會(huì)變差。

傳統(tǒng)的Hessian矩陣近似方法中,當(dāng)構(gòu)建第i行時(shí),只考慮雅可比矩陣第i列的自相關(guān)和第i列與相鄰幾列的互相關(guān)。傳統(tǒng)Hessian近似矩陣的結(jié)構(gòu)如圖3(a)所示。對(duì)角線元素是圖3(b)中矩形區(qū)域中心網(wǎng)格的自相關(guān),副對(duì)角線元素是中心網(wǎng)格與其他相鄰網(wǎng)格的互相關(guān)。對(duì)比圖2(b)和圖3(b),認(rèn)為圖2(b)的近似Hessian矩陣更精確,因?yàn)樗俗钪匾脑?而且還包括了z方向的互相關(guān),在傳統(tǒng)近似方法中該信息被忽略。Ha和Shin[2]在構(gòu)建Hessian近似矩陣時(shí)需要上百個(gè)非零元素,SPNLCG在對(duì)搜索方向步長(zhǎng)進(jìn)行優(yōu)化后,僅需要5個(gè)非零元素。

與NLCG相比SP-NLCG唯一額外的運(yùn)算是求解式(12)左邊第一項(xiàng),令Cigi=,因?yàn)榻艸essian矩陣的高度稀疏性,采用了LU分解直接求解近似Hessian矩陣的逆矩陣。采用稀疏預(yù)條件算子SP-NLCG的效率與NLCG相當(dāng),但是性能卻得到明顯改善。

圖2 稀疏近似Hessian矩陣結(jié)構(gòu)及其構(gòu)建網(wǎng)格示意圖Fig.2 Structure of sparse approximate Hessian matrix and its diagram of inversion grids constructing matrix

圖3 傳統(tǒng)近似Hessian矩陣結(jié)構(gòu)及其構(gòu)建網(wǎng)格示意圖Fig.3 Structure of traditional approximate Hessian matrix and its diagram of inversion grids constructing matrix

2.4 正則化參數(shù)自適應(yīng)更新方法

Rodi等[1]的研究中對(duì)正則化參數(shù)λ沒有進(jìn)行討論,只是將其固定為一個(gè)正實(shí)數(shù),當(dāng)預(yù)條件算子為Ci=(γiI+λLTL)-1時(shí),因?yàn)棣胕是一個(gè)正實(shí)數(shù),LTL是對(duì)稱正定矩陣,所以C總是可逆,λ的作用僅為控制模型的平滑度,λ越大求得的模型越平滑。然而當(dāng)用稀疏近似Hessian矩陣代替預(yù)條件算子的單位矩陣時(shí),γiH不一定可逆,需要λLTL項(xiàng)改善矩陣的條件數(shù),若λ過大則降低了模型的分辨率且增加迭代次數(shù),若λ過小,則γiH+λLTL不可逆,導(dǎo)致計(jì)算失敗。因此,理想的方案是在迭代過程中自適應(yīng)修改λ,既能保持模型的分辨率,又能保證算法的穩(wěn)定性和計(jì)算效率。

本文中采用的自適應(yīng)更新正則化參數(shù)λ的方法如下:選擇正則化因子序列λ0＞λ1＞＞λ2＞…＞＞λl,其中λ0要足夠大而λl要足夠小,在反演過程中根據(jù)每次迭代情況從序列中選擇下次迭代的正則化參數(shù)。正則化參數(shù)序列可以用下式計(jì)算:λk=λ0qk,k=1,2,…,l;q＜0。

正則化參數(shù)自適應(yīng)更新算法的計(jì)算流程如下:

(1)計(jì)算當(dāng)前模型的目標(biāo)函數(shù)中的穩(wěn)定因子項(xiàng)s(mn);

(2)令ζ=s(mn)/s(mn-1);

(3)如果ζ≤1,則λn+1=λn,如果ζ＞1,則λn+1=λn/ζ;

(4)如果當(dāng)前迭代的數(shù)據(jù)殘差rn減小不夠快, 即0.01,更新正則化參數(shù)

λ′n+1=qλn+1,衰減系數(shù)q通常在區(qū)間[0.5,0.9]上選取。

正則化參數(shù)的初始值λ0計(jì)算公式為

式中,mapr為先驗(yàn)?zāi)Ｐ?m0為初始模型;F為正演算子。通過自適應(yīng)算法,實(shí)現(xiàn)了反演精度和計(jì)算穩(wěn)定性之間的平衡。

3 數(shù)值實(shí)驗(yàn)

3.1 塊體組合模型

圖4(a)所示模型與文獻(xiàn)[7]中模型相同,均勻半空間中嵌入兩個(gè)矩形異常體,異常體的電阻率分別為5和2000 Ω·m,背景電阻率為100 Ω·m,用此模型生成模擬數(shù)據(jù)。觀測(cè)點(diǎn)11個(gè),間距10 km,頻點(diǎn)16個(gè),頻率范圍為0.001～100 Hz。

用該簡(jiǎn)單模型對(duì)比分析SP-NLCG法與OCCAM法、NLCG法之間的差異。圖4(b)～(d)是不同方法的反演結(jié)果。使用TE和TM數(shù)據(jù)進(jìn)行聯(lián)合反演,參數(shù)為正則化參數(shù)初值λ0=100 000,SP -NLCG的正則化參數(shù)衰減系數(shù)q=0.7,均方根誤差閾值設(shè)為1.3。反演剖面兩側(cè)與底部的不均勻網(wǎng)格是為了減小邊界效應(yīng)所擴(kuò)充的網(wǎng)格。

圖4 塊體組合模型實(shí)驗(yàn)結(jié)果Fig.4 Block model test results

圖4(b)是OCCAM法反演結(jié)果,OCCAM法實(shí)質(zhì)上是正則化的高斯牛頓法,本實(shí)驗(yàn)中用其反演結(jié)果作為標(biāo)準(zhǔn)與其他方法的反演結(jié)果進(jìn)行對(duì)比。圖4(c)為NLCG反演結(jié)果,預(yù)條件算子為Ci=(γiI+ λLTL)-2,沒有使用Hessian矩陣中的元素作為預(yù)條件算子,并且在迭代過程中不更新正則化參數(shù)λ。反演進(jìn)行了15次迭代后因?yàn)椴荒茉俳档途礁`差而終止,此時(shí)均方根誤差為3.17。這是因?yàn)檎齽t化參數(shù)初值過大,且該值在反演過程中保持不變,導(dǎo)致反演結(jié)果過于平滑,分辨率較低。圖4(d)是采用SP-NLCG的反演結(jié)果,預(yù)條件算子為Ci=(γiH+ λiLTL)-2,Hessian圓半徑為1,其效果與OCCAM法相當(dāng),不同之處在于(d)的模型較(b)的粗糙度稍大一些,這是由于兩種方法對(duì)正則化參數(shù)的更新策略不同所導(dǎo)致的,OCCAM法在每次迭代中為了尋找最優(yōu)的正則化參數(shù)要進(jìn)行多次正演計(jì)算,在本實(shí)驗(yàn)中每次迭代需要6～8次正演計(jì)算尋優(yōu)正則化參數(shù),從而導(dǎo)致計(jì)算量大幅度增加。本文中采用的正則化參數(shù)更新算法不需要額外的正演計(jì)算,盡管SP -NLCG達(dá)到收斂所需的迭代次數(shù)(60次)多于OCCAM法(12次),但耗時(shí)卻比OCCAM法少。

本實(shí)驗(yàn)?zāi)Ｐ头浅：?jiǎn)單,采用不同Hessian圓半徑的反演結(jié)果都非常接近。圖5為不同反演方法迭代誤差曲線,從圖中可以看出OCCAM法收斂速率最快,因此它需要的迭代次數(shù)最少。SP-NLCG最終也收斂到誤差閾值,當(dāng)Hessian圓半徑為0時(shí),即預(yù)條件算子僅取Hessian矩陣對(duì)角線元素,所需迭代次數(shù)大于Hessian圓半徑不為0時(shí)的迭代次數(shù)。當(dāng)Hessian圓半徑為1時(shí),比半徑為3和5時(shí)所需迭代次數(shù)少,計(jì)算速度最快。說明對(duì)于SP-NLCG,預(yù)條件算子中包含Hessian矩陣元素并非越多越好,本實(shí)驗(yàn)以及后續(xù)實(shí)驗(yàn)證明當(dāng)Hessian圓半徑為1,即預(yù)條件算子每行包含5個(gè)Hessian矩陣元素時(shí),反演性能最優(yōu)。

圖5 塊體組合模型反演誤差曲線Fig.5 RMS curves of block model test

3.2 推覆逆掩模型實(shí)驗(yàn)

為了測(cè)試SP-NLCG法對(duì)復(fù)雜構(gòu)造的適用性,設(shè)計(jì)一個(gè)推覆逆掩模型(圖6(a))。該模型的右側(cè)為沉積地層,電阻率從淺到深分別為100、30、100、400 Ω· m的老地層逆掩在新地層之上。觀測(cè)點(diǎn)34個(gè),間距1 km,頻點(diǎn)數(shù)為38,頻率范圍為0.001～320 Hz。

圖6(b)～(e)是不同方法的反演結(jié)果,使用TE和TM數(shù)據(jù)進(jìn)行聯(lián)合反演,正則化參數(shù)初值λ0= 100000,SP-NLCG的正則化參數(shù)衰減系數(shù)q=0.7,均方根誤差閾值設(shè)為1.3。

圖6(c)為NLCG反演結(jié)果,反演進(jìn)行了14次迭代后因?yàn)椴荒茉俳档途礁`差而終止,此時(shí)均方跟誤差為5.69。與圖6(b)的OCCAM法反演結(jié)果相比,反演精度明顯要低很多,與圖6(d)所示的Hessian圓半徑為1的SP-NLCG的15次迭代結(jié)果相比,其精度也低很多。原因有兩方面:①預(yù)條件算子中不包含任何Hessian矩陣元素,損失大量有效信息,降低了反演分辨率;②反演迭代過程中不更新正則化參數(shù),過大的正則化參數(shù)初值使得模型過于平滑,后期迭代已經(jīng)不能再降低均方根誤差,導(dǎo)致反演不能達(dá)到規(guī)定閾值而提前結(jié)束。

圖6(e)為Hessian圓半徑為1的SP-NLCG的36次迭代結(jié)果,其精度與OCCAM法反演結(jié)果相當(dāng),很好地恢復(fù)出了逆掩構(gòu)造和沉積地層。Hessian半徑取不同值時(shí)的反演結(jié)果都比較接近,但計(jì)算效率存在較大差異,圖7給出了不同反演方法迭代誤差曲線。

圖6 推覆逆掩模型實(shí)驗(yàn)結(jié)果Fig.6 Overthrust model test results

從圖7和表1中可以看出SP-NLCG的收斂速度介于OCCAM法和NLCG法之間,當(dāng)Hessian圓半徑為0時(shí),反演所需45次迭代達(dá)到誤差閾值, Hessian圓半徑為5時(shí),需要40次迭代,Hessian圓半徑為1和3時(shí),分別需要36次和35次迭代,盡管Hessian圓半徑為1時(shí)的迭代次數(shù)比Hessian圓半徑為3時(shí)迭代次數(shù)要多一次,然而從表1中可以看到其計(jì)算時(shí)間卻是最短的,這是因?yàn)轭A(yù)條件算子包含的非零元素越多,對(duì)其進(jìn)行求逆運(yùn)算所需的時(shí)間越長(zhǎng)。

圖7 推覆逆掩模型反演誤差曲線Fig.7 RMS curves of overthrust model test

圖8是SP-NLCG法在Hessian圓選取不同半徑時(shí)預(yù)條件算子包含的Hessian矩陣元素,其中圖8(a)是完整的Hessian矩陣,圖8(b)～(d)分別是Hessian圓半徑為0、1、3時(shí)的近似Hessian矩陣。從圖中可以看出近似Hessian矩陣都為高度稀疏的對(duì)角陣,元素集中在主副對(duì)角線上,元素個(gè)數(shù)比完整Hessian矩陣少很多,因此在求預(yù)條件算子時(shí)無論是計(jì)算量還是對(duì)內(nèi)存的要求都不高。

表1 不同反演方法性能參數(shù)對(duì)比Table 1 Performance parameters of different inversion methods

圖8 Hessian矩陣對(duì)比Fig.8 Hessian matrix comparison

4 實(shí)測(cè)資料處理

為了檢驗(yàn)本文研究方法處理實(shí)測(cè)數(shù)據(jù)的有效性,對(duì)中國(guó)西部某地的實(shí)測(cè)MT數(shù)據(jù)進(jìn)行反演。該地區(qū)屬于山前構(gòu)造帶,第四系為不穩(wěn)定中高電阻層;第三系到白堊系為低阻層,電阻率為幾歐姆米;侏羅系為次低阻層,電阻率一般小于30 Ω·m;三疊系為次高阻層,電阻率在30～100 Ω·m范圍變化;石炭系則基本是高阻,超過100歐姆甚至達(dá)上千歐姆米,也有部分上石炭統(tǒng)電阻率不到100 Ω·m。

圖9上圖為MT 2維連續(xù)介質(zhì)反演剖面,剖面總長(zhǎng)60 km,觀測(cè)點(diǎn)120個(gè),觀測(cè)點(diǎn)距500 m,觀測(cè)數(shù)據(jù)頻點(diǎn)數(shù)為38,頻率范圍為0.001～320 Hz。剖面能較清晰地反映沉積盆地的基本構(gòu)造形態(tài),但是對(duì)深部構(gòu)造的恢復(fù)不理想,出現(xiàn)了許多不自然的“掛面條”現(xiàn)象,難以刻畫深層的盆山接觸關(guān)系。剖面中下部偏左處的相對(duì)低阻隱約反映出存在逆掩構(gòu)造,但是該構(gòu)造在地震偏移剖面(圖9下圖)很難識(shí)別,因此很難確定該構(gòu)造的存在。

圖10為SP-NLCG法反演得到的電阻率剖面和電震聯(lián)合解釋剖面。可以看出,該剖面對(duì)沉積地層的形態(tài)刻畫清晰,且恢復(fù)出了深層的逆掩構(gòu)造。最上面的低阻地層為大面積出露地表的第四系-白堊系,從左到右逐漸變薄,最右邊厚度為1800 m;第二層為侏羅系-三疊系的相對(duì)中高阻層,被石炭系高阻推覆體沖斷成東西兩段,厚度約2 600 m;第三層和第四層分別為上二疊統(tǒng)的中低阻層和下二疊統(tǒng)的中高阻層,這兩層都被逆沖斷層劃成3段,第五層為石炭系的高阻層,在測(cè)線最右邊有出露,剖面中部為大套的高阻層,并未發(fā)育低阻層,根據(jù)地震偏移剖面解釋為石炭和二疊地層。利用本文方法的反演結(jié)果指導(dǎo)了該地區(qū)構(gòu)造解釋,尤其是對(duì)深層部位石炭系和二疊系的刻畫較為精準(zhǔn),證實(shí)了SP-NLCG法解決中國(guó)西部復(fù)雜構(gòu)造問題的有效性。

圖9 MT 2維連續(xù)介質(zhì)反演剖面與地震偏移剖面Fig.9 MT 2D continuous inversion result and seismic migration profile

圖10 預(yù)條件非線性共軛梯度法反演剖面與電震聯(lián)合解釋剖面Fig.10 Inversion results of SP-NLCG and MT&seismic co-interpretive profile

5 結(jié) 論

(1)通過分析雅克比矩陣中列向量的自相關(guān)和互相關(guān),發(fā)現(xiàn)構(gòu)建近似Hessian矩陣最有效的模式是僅使用以某網(wǎng)格為中心的Hessian圓內(nèi)的網(wǎng)格所對(duì)應(yīng)的雅克比矩陣列向量。數(shù)值實(shí)驗(yàn)證明該方法的反演結(jié)果精度比非線性共軛梯度法精度更高,與高斯牛頓法相當(dāng),且計(jì)算速度比高斯牛頓法快,與非線性共軛梯度法相當(dāng),在通常情況下Hessian圓半徑為1時(shí)反演性能最優(yōu)化。

(2)通過引入正則化參數(shù)的自適應(yīng)更新算法,在迭代過程中自適應(yīng)修改λ,既能保持模型的分辨率,又能保證算法的穩(wěn)定性和計(jì)算效率。

[1] RODI W,MACKIE R.Nonlinear conjugate gradients algorithm for 2D magnetotelluric inversion[J].Geophysics,2001,66:174-187.

[2] HA T,SHIN C.Problems of constructing the approximate Hessian in waveform inversion[J].SEG Technical Program Expanded Abstracts,2002,21:961-962.

[3] MACKIE R L,BENNET B R,MADDEN T R.Long-period magnetotelluric measurements near the central California coast:a land-locked view of conductivity structure under the Pacific Ocean[J].Geophysical Journal International,1988,95:181-194.

[4] 金建銘,王建國(guó).電磁場(chǎng)有限元方法[M].西安:西安電子科技大學(xué)出版社,1998.

[5] TIKHONOV.Solution of incorrectly formulated problems and the regularization method[J].Soviet Math,1963,4: 1035-1038.

[6] TIKHONOV.Regularization of incorrectly posed problems [J].Soviet Math,1963,4:1624-1627.

[7] CONSTABLE S C,PARKER R L,CONSTABLE C G. Occam's inversion:a practical algorithm for generating smooth models from electromagnetic sounding data[J]. Geophysics,1987,52:289-300.

[8] DEGROOT HEDLIN C,CONSTABLE S.Occam's inversion to generate smooth,two-dimensional models from magnetotelluric data[J].Geophysics,1990,55(12): 1613-1624.

[9] PORTNIAGUINE O,ZHDANOV M S.Focusing geophysical inversion images[J].Geophysics,1999,64:874-887.

[10] PORTNIAGUINE O,ZHDANOV M S.3-D magnetic inversion with data compression and image focusing[J]. Geophysics,2002,67(5):1532-1541.

[11] ZHDANOV M S,ELLIS R,MUKHERJEE S.Three-D regularized focusing inversion of gravity gradient tensor data[J].Geophysics,2004,69:925-937.

[12] FLETCHER R,REEVES C M.Function minimization by conjugate gradients[J].Computer J,1959,7:149-154.

(編輯 修榮榮)

Two-dimensional MT inversion method based on an improved preconditioned nonlinear conjugate gradient algorithm

XIANG Peng
(Research Institute of New District in West China,Shengli Oilfield,Dongying 257000,China)

Accuracy and computational efficiency are conflicting to each other in MT inversion methods.Gauss-Newton methods have high accuracy but poor efficiency,while nonlinear conjugate gradient methods are computationally efficient but less accurate.A new preconditioned nonlinear conjugate gradient method is proposed to improve both accuracy and computational cost by using a precondition operator whose performance is approximated by Gauss-Newton Hessian matrix.A self-adapted algorithm updating the regularization parameter is used to guarantee the balance between the stability and the accuracy of the inversion.Through the model tests and comparisons with other methods,the proposed method is proved achieving the goal. The method is applied to invert the MT data measured from a region in west China.The inverted resistivity profile is conjointly interpreted with seismic migration profile,and illustrates its significance in complicated structural interpretation.

magnetotelluric(MT);preconditioned nonlinear conjugate gradient(PNLCG);regularization parameter;self-adapted algorithm;inversion

P 631.3

A

1673-5005(2014)04-0042-08

10.3969/j.issn.1673-5005.2014.04.006

2013-12-19

國(guó)家科技重大專項(xiàng)(2011ZX05002-002)

相鵬(1978-),男,工程師,博士,主要從事地球物理正反演方法研究。E-mail:ddzk@163.com。

相鵬.一種改進(jìn)的二維MT預(yù)條件非線性共軛梯度反演方法[J].中國(guó)石油大學(xué)學(xué)報(bào):自然科學(xué)版,2014, 38(4):42-49.

XIANG Peng.Two-dimensional MT inversion method based on an improved preconditioned nonlinear conjugate gradient algorithm[J].Journal of China University of Petroleum(Edition of Natural Science),2014,38(4):42-49.