一種利用特征值性質(zhì)的MT阻尼最小二乘反演

唐榮江,王緒本,2,甘 露

(1.成都理工大學(xué)地球物理學(xué)院,四川成都610059;2.地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,四川成都610059)

唐榮江,王緒本,甘露.一種利用特征值性質(zhì)的MT阻尼最小二乘反演[J].石油物探,2017,56(6):-904

TANG Rongjiang,WANG Xuben,GAN Lu.A damped least square inversion for MT utilizing eigenvalue property[J].Geophysical Prospecting for Petroleum,2017,56(6):-904

一種利用特征值性質(zhì)的MT阻尼最小二乘反演

唐榮江1,王緒本1,2,甘 露1

(1.成都理工大學(xué)地球物理學(xué)院,四川成都610059;2.地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,四川成都610059)

阻尼最小二乘反演的方法總會(huì)得到一個(gè)最終迭代方程,該方程系數(shù)矩陣的特征值對(duì)反演結(jié)果有重要影響:小特征值相對(duì)大特征值對(duì)模型修正量的貢獻(xiàn)更大,能夠得到更高的分辨率,但同時(shí)更加不穩(wěn)定。為此,以一維層狀介質(zhì)的大地電磁阻尼最小二乘反演為基礎(chǔ),采用特定單位的電阻率參數(shù)和層厚參數(shù)進(jìn)行運(yùn)算,同時(shí)改進(jìn)阻尼因子的自適應(yīng)調(diào)節(jié)方法,可以巧妙地利用小特征值能提高模型空間分辨率的特點(diǎn),得到一種改進(jìn)的阻尼最小二乘反演方法,并與傳統(tǒng)阻尼最小二乘法進(jìn)行了對(duì)比。合成數(shù)據(jù)測(cè)試結(jié)果表明,該方法能夠得到電性陡變的反演模型,彌補(bǔ)了平滑模型反演不能有效識(shí)別電性突變界面的不足,對(duì)電性變化較為劇烈的地電結(jié)構(gòu)的識(shí)別具有較好效果,分辨率有一定提高。

阻尼最小二乘反演;大地電磁;特征值;阻尼因子;分辨率;突變界面

阻尼最小二乘法最早由LEVENBERG提出,1963年MARQUARDT[1]對(duì)其做了完善,1971年FLETCHER[2]做了進(jìn)一步改進(jìn),改進(jìn)的阻尼最小二乘法是目前求解無(wú)約束最優(yōu)化問(wèn)題最優(yōu)秀的算法之一。1988年陳鐘琦[3]討論了高阻尼因子的影響并提出了改進(jìn)方案。傳統(tǒng)的大地電磁反演方法大都是基于最大平滑理論,反演出的地電圖像都是漸變的,不能準(zhǔn)確反映出電性突變界面的位置。例如CONSTABLE等[4]提出的Occam反演,RODI等[5]提出的非線性共軛梯度反演(nonlinear conjugate gradient inversion,NLCG)都是光滑模型反演。

最小二乘反演將極小值問(wèn)題轉(zhuǎn)化成求解線性方程組,MARQUARDT[1]提出阻尼最小二乘法并改進(jìn)了方程組,增大了矩陣的特征值,降低了條件數(shù),改善了收斂性。奇異值分解法仍是目前求病態(tài)解線性方程組最有力的工具之一,并在近幾年發(fā)展了基于奇異值分解法的地震數(shù)據(jù)隨機(jī)噪聲壓制方法和子波提取方法[6-7]。根據(jù)奇異值分解理論,WIGGINS[8]曾用截?cái)嘁恍┬∑娈愔档姆椒▉?lái)減少解,以取得與阻尼最小二乘法類似的效果,但是應(yīng)當(dāng)截除什么樣的奇異值是奇異值截?cái)喾ㄖ械闹饕щy,且截?cái)嘈∑娈愔岛蠓直媛室搽S之降低。李平等[9]探討了奇異值分解中廣義逆解的可靠性估計(jì),給出了更加合理的奇異值修改方案。

阻尼最小二乘法最終迭代矩陣為滿秩對(duì)稱方陣,相對(duì)于雅可比矩陣來(lái)說(shuō)形式更簡(jiǎn)易,且雅可比矩陣的奇異值和迭代方陣的特征值對(duì)反演結(jié)果具有相同影響。阻尼因子的大小直接影響矩陣特征值的大小,小特征值與高分辨率對(duì)應(yīng),大特征值與較小的方差和較好的收斂性對(duì)應(yīng)。如果能有效利用和改進(jìn)小特征值,同時(shí)注意大特征值以保證方程的收斂性,是否能在一定程度上提高阻尼最小二乘法的分辨率,或針對(duì)不同的地質(zhì)模型采取不同的特征值變化方案以達(dá)到約束反演的目的？

通常情況下,大地電磁一維反演的分辨率要高于二維、三維反演的分辨率,馬奎特一維反演中層厚和電阻率要同時(shí)參與運(yùn)算,這為提高一維反演的分辨率和優(yōu)化反演結(jié)果提供了更大的改進(jìn)空間。據(jù)此,本文以一維層狀模型的大地電磁阻尼最小二乘反演為基礎(chǔ),進(jìn)行特征值分析,探討其如何影響模型空間,據(jù)此提出改進(jìn)方法:不直接修改特征值,而是使用特定單位的反演參數(shù)和改變阻尼因子的自適應(yīng)方案以間接達(dá)到改變特征值的目的,從而優(yōu)化反演結(jié)果。理論數(shù)據(jù)實(shí)驗(yàn)表明,改進(jìn)后的反演方法能夠有效識(shí)別電性突變界面,彌補(bǔ)了光滑模型反演的缺陷,一定程度上提高了大地電磁的分辨率。

1 反演基本原理

大地電磁阻尼最小二乘法的目標(biāo)函數(shù)為:

(1)

式中:d為觀測(cè)數(shù)據(jù);α為阻尼因子;m=[ρ1…ρj…ρN,h1…h(huán)j…h(huán)N-1]為模型參數(shù),其中,ρj為第j層電阻率,hj為第j層厚度,N為反演層數(shù);A(m)為正演算子;Wd和Wm分別為數(shù)據(jù)空間和模型空間加權(quán)矩陣。為了將非線性問(wèn)題線性化,考慮:

(2)

用泰勒公式對(duì)A(m1)進(jìn)行一階展開(kāi),F為A的弗雷歇導(dǎo)數(shù),若Wd和Wm為單位矩陣,A(m1)表示為A(m0)+Fm0Δm,目標(biāo)函數(shù)變?yōu)?

(3)

求目標(biāo)函數(shù)對(duì)ΔmT的導(dǎo)數(shù):

(4)

根據(jù)最小化條件要求導(dǎo)數(shù)為0,可得:

(5)

通過(guò)解方程(5),可得到Δm,多次迭代之后可得到滿足精度的反演解。在大地電磁反演中,為了增強(qiáng)方程的穩(wěn)定性,陳樂(lè)壽等[10]對(duì)目標(biāo)函數(shù)做了改進(jìn):對(duì)數(shù)據(jù)空間取對(duì)數(shù),并對(duì)模型空間取比例。據(jù)此,方程(1)變?yōu)?

(6)

(7)

在常規(guī)的阻尼最小二乘法中,為了避免電阻率本身作為模型參數(shù)數(shù)量級(jí)差異太大造成反演不穩(wěn)定,反演對(duì)電阻率和層厚作為模型參數(shù)時(shí)取其對(duì)數(shù):

(8)

那么雅可比矩陣中:

根據(jù)上述理論,常規(guī)阻尼最小二乘反演具體實(shí)現(xiàn)步驟[10]如下。

1) 設(shè)置初始模型m,初始阻尼因子α(這里設(shè)置為1),加權(quán)系數(shù)k(通常設(shè)置為10),給定數(shù)據(jù)擬合誤差的閾值E。

4) 判斷條件Φ1(m+Δm)<Φ1(m)是否滿足,如果是,m=m+Δm,α=αk-1,執(zhí)行步驟5),否則α=αk,返回步驟3)。

5) 判斷條件‖lnA(m+Δm)-lnd‖

2 特征值的利用與反演方法的改進(jìn)

2.1 特征值分析

反演迭代方程的特征值對(duì)反演結(jié)果具有重要影響,下面分析特征值的大小如何影響模型參數(shù)。

(11)

(12)

Λ的具體表達(dá)式為:

(13)

式中:λi為H+αD的特征值。根據(jù)(7)式,模型修正量Δm的求和形式為:

(14)

(15)

式中:N為層厚;M為測(cè)點(diǎn)數(shù)。由(14)式可知,特征值λi越小,對(duì)模型參數(shù)Δmj的影響越大,同時(shí)也使方程更加不穩(wěn)定。

2.2 反演參數(shù)的改進(jìn)

反演參數(shù)的改進(jìn)為:避免取對(duì)數(shù),直接使用以千米為單位的層厚和歐姆米為單位的電阻率進(jìn)行反演。這種改進(jìn)巧妙地利用了小特征值對(duì)模型修正量貢獻(xiàn)更大的特點(diǎn),使得層厚修正量比例大于電阻率修正量比例,以達(dá)到持續(xù)改變反演模型的目的。

2.3 阻尼因子的改進(jìn)與反演的收斂性

實(shí)際計(jì)算中,阻尼因子α不能取太大,否則步長(zhǎng)太小,收斂太慢。僅當(dāng)Φ1(m+Δm)>Φ1(m)時(shí),才被迫取較大的阻尼因子,使迭代穩(wěn)定收斂。當(dāng)α(0)給定后,逐步嘗試增加或減少。若采用本文方法的反演參數(shù),層厚數(shù)值較小,在α較小時(shí)方程穩(wěn)定性降低,層厚參數(shù)極易變?yōu)樨?fù)值,使反演失敗。此外,根據(jù)2.2節(jié),為了使Δρ與Δh盡快接近以達(dá)到層厚改變的目的,必須保證阻尼因子的快速下降?；谝陨夏康?我們改進(jìn)阻尼因子的自適應(yīng)變化方案,反演的步驟4)變?yōu)?

4)m=m+Δm,判斷各層的mi是否大于0,如果是,則α=αk-1,返回步驟2),否則α=αk,返回步驟3)。

由表1可知:①本文方法的α在第7次迭代之后就遠(yuǎn)遠(yuǎn)低于傳統(tǒng)方法的α,說(shuō)明新的收斂條件可以保證α更快的下降速度;②兩種方法的擬合差均逐步減少,本文方法的擬合差沒(méi)有傳統(tǒng)方法的小,但阻尼因子更低,導(dǎo)致特征值更低,觀測(cè)數(shù)據(jù)的微小變化將導(dǎo)致模型參數(shù)的更大變化;③本文方法中第2次迭代到第3次迭代擬合差不降反升,且后期震蕩變化,整體為降低的趨勢(shì),說(shuō)明阻尼因子起著調(diào)節(jié)特征值和穩(wěn)定

表1 不同迭代次數(shù)下阻尼因子和擬合差變化規(guī)律

性的作用,當(dāng)層厚修正量過(guò)大使層厚出現(xiàn)負(fù)值時(shí),則提高阻尼因子,降低層厚修正量,保證方程穩(wěn)定性;若層厚為正值,則降低阻尼因子,保證分辨率。以上結(jié)論說(shuō)明通過(guò)判斷模型參數(shù)是否為負(fù)值來(lái)保證反演的穩(wěn)定性和收斂性是合理的。

2.4 分辨率提高的原因

前文的理論說(shuō)明:采用特定的反演參數(shù)單位保證了反演迭代方程具有更小的特征值,使得層厚參數(shù)相對(duì)于電阻率有更高的修正量;本文方法保證了反演的穩(wěn)定性和收斂性,要說(shuō)明的是,新的改進(jìn)能夠達(dá)到更高的分辨率,迭代后期層厚參數(shù)修正量必須以負(fù)值為主,使平滑的反演模型得到壓縮,突出異常。但使用定量分析證明模型修正量為負(fù)值較為困難,可作定性分析:在保證數(shù)據(jù)擬合差極小的情況下,相對(duì)于增加層厚,減少層厚來(lái)保證數(shù)據(jù)的高度擬合明顯要容易得多。為了說(shuō)明其正確性,這里同樣采用2.3節(jié)中的模

型和參數(shù),使用本文改進(jìn)方法對(duì)不同迭代次數(shù)的模型修正量成圖,并對(duì)比兩種方法的反演結(jié)果。

在圖1中,圖1a和圖1b的縱坐標(biāo)分別為模型修正量Δm和模型修正量Δm與模型參數(shù)m(包括電阻率ρ和層厚h)之比,橫坐標(biāo)為反演參數(shù)序號(hào),圖1a和圖1b 中的反演參數(shù)序號(hào)1～15對(duì)應(yīng)第1層到15層的電阻率,16～29對(duì)應(yīng)第1層到14層的層厚。由圖可見(jiàn),圖1a中Δh遠(yuǎn)小于Δρ,且兩者在迭代后期都趨于0;圖1b中Δh/h在早期總體小于Δρ/ρ,隨后逐漸增大甚至大于Δρ/ρ,兩者在后期都減小,但Δh/h總體仍大于Δρ/ρ,且主要為負(fù)值。應(yīng)當(dāng)注意,引起層厚變化率大于電阻率變化率的并不是模型修正量本身,而是修正量與模型參數(shù)的比例。

圖2a和圖2b分別為對(duì)單個(gè)低阻層模型采用傳統(tǒng)方法和本文方法進(jìn)行20次迭代時(shí)的反演結(jié)果,可以看出,相對(duì)于本文方法,傳統(tǒng)方法的反演結(jié)果較為平緩,不能準(zhǔn)確反映電性突變界面的地下結(jié)構(gòu)。從圖2b 可以看出,電性過(guò)度界面節(jié)點(diǎn)較為密集,說(shuō)明該處的層厚受到壓縮,在一定程度上提高了分辨率,同時(shí)表明,利用層厚的小特征值可以提高模型空間的分辨率,而傳統(tǒng)方法擬合差更低。該結(jié)論與小特征值對(duì)應(yīng)高分辨率,大特征值意味著以犧牲一些分辨率為代價(jià)來(lái)?yè)Q取較小的方差的結(jié)論一致。

圖1 不同迭代次數(shù)下采用本文反演方法的模型修正量 a 模型修正量; b 模型修正量與模型參數(shù)比值

圖2 單個(gè)低阻層模型傳統(tǒng)方法(a)與本文方法(b)進(jìn)行20次迭代的反演結(jié)果

本文方法通過(guò)壓縮電性界面突變區(qū)層厚來(lái)突出異常,所以對(duì)于平滑模型的反演,本文方法得出的結(jié)果往往與實(shí)際情況不吻合,此時(shí)選用傳統(tǒng)阻尼最小二乘反演能取得更為理想的結(jié)果。所以,在實(shí)際情況中,針對(duì)不同的反演目標(biāo),應(yīng)當(dāng)選擇合適的反演方法來(lái)解決地質(zhì)問(wèn)題。

2.5 初始模型和解方程組

從本質(zhì)上來(lái)說(shuō),本文方法仍然是阻尼最小二乘法的延續(xù),所以本文方法和傳統(tǒng)阻尼最小二乘反演法的初始模型、方程組解法、收斂條件等的選擇原則一致。本文解方程組采用共軛梯度法[12],該方法可以快速、穩(wěn)定求解對(duì)稱矩陣的線性代數(shù)方程組;收斂條件上可以采用設(shè)定閾值或者固定迭代次數(shù)的方式,反演深度和加權(quán)系數(shù)在各理論模型的反演中均有說(shuō)明。

初始模型建立可以采用均勻半空間和Bostic變換,本文采用類似于均勻半空間的改進(jìn)方法:以大地電磁最高頻視電阻率當(dāng)作初始模型第1層電阻率,將最低頻視電阻率當(dāng)作初始模型最后一層電阻率,其它層的電阻率位于兩者之間漸變分布,各層厚度采用等比增長(zhǎng)。

3 合成數(shù)據(jù)測(cè)試及分析

設(shè)計(jì)了2個(gè)理論模型來(lái)驗(yàn)證本文算法的有效性。圖3的6層理論模型與文獻(xiàn)[13]中的模型一致,圖4的8層理論模型與文獻(xiàn)[14]中的模型一致,無(wú)隨機(jī)噪聲,反演模型設(shè)定為35層,6層模型反演深度15km,8層模型反演深度為40km,阻尼因子初始值都為10,加權(quán)系數(shù)k為10,迭代次數(shù)為20次。

圖3 采用本文方法和傳統(tǒng)阻尼最小二乘法的6層模型反演結(jié)果對(duì)比 a 傳統(tǒng)方法; b 本文方法; c 擬合曲線

圖4 采用本文方法和傳統(tǒng)阻尼最小二乘法的8層模型反演結(jié)果對(duì)比 a 傳統(tǒng)方法; b 本文方法; c 擬合曲線

圖3和圖4分別顯示了采用本文方法、傳統(tǒng)阻尼最小二乘法的6層模型和8層模型的反演結(jié)果。由圖3和圖4可見(jiàn),采用傳統(tǒng)方法能基本反演出各層地電模型,吻合程度相對(duì)較低,電性突變界面的反演層過(guò)渡相對(duì)平緩,不能準(zhǔn)確識(shí)別突變界面的具體位置(圖3a,圖4a);采用本文方法不但能有效反演出各層地電模型,而且在很大程度上與理論模型吻合,能有效識(shí)別電性突變界面,相對(duì)傳統(tǒng)方法來(lái)說(shuō),分辨率有一定提高(圖3b,圖4b);兩種算法都能較好地?cái)M合原始數(shù)據(jù),雖然傳統(tǒng)方法具有更高的擬合精度,但從圖像上無(wú)法直接看出,這說(shuō)明可以接受本文算法的擬合精度,且一定程度上犧牲了擬合精度來(lái)提高分辨率是值得的(圖3c,圖4c)。

4 結(jié)束語(yǔ)

本文以迭代方程特征值的性質(zhì)為理論基礎(chǔ),進(jìn)行反演參數(shù)和阻尼因子自適應(yīng)方法的相關(guān)改進(jìn),并將改進(jìn)后的大地電磁一維阻尼最小二乘算法有效地應(yīng)用于解決電性突變地質(zhì)結(jié)構(gòu)的反演問(wèn)題。研究表明,小特征值對(duì)應(yīng)大的模型修正量和更高的分辨率,大特征值對(duì)應(yīng)小的模型修正量和更小的擬合差;較大的特征值差異導(dǎo)致較大模型修正量差異,較小的特征值差異導(dǎo)致較小模型修正量差異。使用了數(shù)量級(jí)懸殊的電阻率(Ω·m)和層厚(km)直接反演,可以導(dǎo)致對(duì)應(yīng)的特征值差異。阻尼因子較大時(shí),Δρ較大,Δh較小,隨著α的降低,Δρ和Δh逐漸接近,數(shù)量級(jí)較小的層厚更易發(fā)生改變,且修正量主要為負(fù)值。同時(shí),采取判斷模型參數(shù)是否全為正來(lái)修正阻尼因子α,力求在保證反演收斂情況下,最快速度下降阻尼因子,達(dá)到使Δρ和Δh逐漸接近的目的。

一維層狀理論模型的反演結(jié)果表明:改進(jìn)后的反演方法能夠有效識(shí)別電性突變界面,針對(duì)地下電性突變的結(jié)構(gòu),在一定程度上提高了大地電磁測(cè)深的分辨率,且該方法可以拓展應(yīng)用到可控源電磁法的反演上,但是對(duì)于地下電性平滑的界面,采用傳統(tǒng)方法反演則更有效。

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

[2] FLETCHER R.A modified Marquardt subroutine for non-linear least squares:UKAEA Research Group Report AERF R-6799[R].UK:Atomic Energy Research Establishment,1971

[3] 陳鐘琦.高阻尼因子對(duì)阻尼最小二乘法效果的影響和克服[J].現(xiàn)代地質(zhì),1988,2(2):255-264

CHEN Z Q.The influence of the high damping factor to the result of damped least square method and it’s overcome[J].Geoscience,1988,2(2):255-264

[4] 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(3):289-300

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

[6] 馬繼濤,王建花,劉國(guó)昌.基于頻率域奇異值分解的地震數(shù)據(jù)插值去噪方法研究[J].石油物探,2016,55(2):205-213

MA J T,WANG J H,LIU G C.Seismic data noise attenuation and interpolation using singular value decomposition in frequency domain[J].Geophysical Prospecting for Petroleum,2016,55(2):205-213

[7] 王蓉蓉,戴永壽,李闖,等.基于奇異值分解的時(shí)變子波提取準(zhǔn)確性評(píng)價(jià)方法[J].石油物探,2015,54(5):531-540

WANG R R,DAI Y S,LI C,et al.An evaluation criterion on the accuracy of time-varying wavelet extraction based on singular value decomposition[J].Geophysical Prospecting for Petroleum,2015,54(5):531-540

[8] WIGGINS R A.The general linear inverse problem:implication of surface waves and free oscillations for earth structure[J].Reviews of Geophysics,1972,10(1):251-285

[9] 李平,王椿鏞,徐厚澤,等.地球物理反演中奇異值分解應(yīng)用的若干問(wèn)題探討[J].自然科學(xué)進(jìn)展,2001,11(8):891-896

LI P,WANG C Y,XU H Z,et al.Several problems discussion in the application of singular value decomposition in geophysical inversion[J].Progress in Natural Science,2001,11(8):891-896

[10] 陳樂(lè)壽,劉任,壬天生.大地電磁測(cè)深資料處理與解釋[M].北京:石油工業(yè)出版社,1989:1-236

CHEN L S,LIU R,REN T S,et al.Magnetotelluric sounding data processing and interpretation[M].Beijing:Petroleum Industry Publishing House,1989:1-236

[11] 廖安平,劉建州.矩陣論[M].長(zhǎng)沙:湖南大學(xué)出版社,2005:97-101

LIAO A P,LIU J Z.Theory of matrices[M].Changsha:Hunan University Publishing House,2005:97-101

[12] NOCEDAL J,WRIGHT S.Numerical optimization 2nd [M].The United States:Springer Science & Business Media,2006:101-133

[13] 管貽亮,李予國(guó),胡祖志,等.大地電磁非線性共軛梯度一維反演[J].石油物探,2014,53(6):752-759

GUAN Y L,LI Y G,HU Z Z,et al.Nonlinera conjugate gradients algorithm for 1D magnetotelluric inversion[J].Geophysical Prospecting for Petroleum,2014,53(6):752-759

[14] 陳小斌,趙國(guó)澤,湯吉,等.大地電磁自適應(yīng)正則化反演算法[J].地球物理學(xué)報(bào),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

附錄A

以下所有矩陣的特征值滿足降序排列λ1>λ2>…>λn。

定理1:若A,B為n階Hermite矩陣,A,B的特征值分別為λi(A),λi(B),則對(duì)于i=1,2,…,n,有:

(A1)

詳細(xì)證明過(guò)程請(qǐng)參見(jiàn)文獻(xiàn)[11],根據(jù)定理1,這里得出如下推論。

推論:若A,B為2N-1階實(shí)對(duì)稱矩陣,B滿足前N個(gè)特征值遠(yuǎn)大于后N-1個(gè)特征值,且遠(yuǎn)大于A的特征值,則有:λi(A+B)的前N個(gè)特征值λs1(B)遠(yuǎn)大于后N-1個(gè)特征值λs2(B);且隨著B(niǎo)矩陣所有特征值成比例逐漸減小,λi(A+B)的前N個(gè)特征值與后N-1個(gè)特征值的差逐漸變小。

證明:

當(dāng)i∈[1,N]時(shí),根據(jù)定理1,不等式左端最大值條件為2N≤r+s≤3N-1,λr(A)+λs(B)要取得最大值,B的特征值必在前N個(gè)中取值;不等式右端最小值條件為1≤r+s≤N+1,λr(A)+λs(B)要取得最小值,由于r≥1,必有s≤N,則B的特征值只能在前N個(gè)中取值。

當(dāng)i∈[N+1,2N-1]時(shí),根據(jù)定理1,不等式左端最大值條件為3N≤r+s≤4N-2,λr(A)+λs(B)要取得最大值,由于r≤2N-1,必有s≥N+1,則B的特征值只能在后N-1個(gè)中取值;不等式右端最小值條件為N+2≤r+s≤2N,λr(A)+λs(B)要取得最小值,由于r≥1,必有s≥N+1,則B的特征值只能在后N-1個(gè)中取值。

綜上可知,λi(A+B)的前N個(gè)特征值的上、下界均包含λs1(B),且1≤s1≤N;λi(A+B)的后N-1個(gè)特征值的上、下界均包含λs2(B),不含λs1(B),且N+1≤s2≤2N-1;所以λi(A+B)的前N個(gè)特征值遠(yuǎn)大于后N-1個(gè)特征值。隨著B(niǎo)矩陣所有特征值成比例逐漸減小,λi(A+B)的前N個(gè)特征值的上、下界大幅減小,因?yàn)棣藄(B)在λr(A)+λs(B)中占有主要比重;后N-1個(gè)特征值的上、下界小幅度降低,因?yàn)闇p小幅度受λr(A)和λs(B)的比例共同控制,所以λi(A+B)的前N個(gè)特征值與后N-1個(gè)特征值的差逐漸變小。

(編輯:顧石慶)

AdampedleastsquareinversionforMTutilizingeigenvalueproperty

TANG Rongjiang1,WANG Xuben1,2,GAN Lu1

(1.CollegeofGeophysics,ChengduUniversityofTechnology,Chengdu610059,China;2.KeyLaboratoryofEarthexplorationandInformationTechniquesofMinistryofEducation,Chengdu610059,China)

Damping least squares inversion always involves a final iterative equation,and the eigenvalues of its coefficient matrix have an important effect on the inversion results.Namely,compared to large eigenvalues,small eigenvalues contribute more to model corrections,resulting in higher resolution,but leading to lower stability.An improved damped least square inversion method is proposed that is based on a one-dimensional magnetotelluric inversion scheme.The method cleverly uses the characteristic that small eigenvalues can improve the spatial resolution,which uses the resistivity and layer thickness parameters of a specific unit for inversion,while improving the adaptive adjustment method of the damping factor.Synthetic data test results showed that,the proposed method could invert a model with an electrical abrupt change,which makes up for the shortcoming of smooth model inversion,which cannot effectively identify an electrical abrupt interface.Compared with traditional methods,this method is more effective for the identification of electrical structures with geoelectrical abrupt change,enabling higher resolution.

damped least square inversion,magnetotelluric (MT),eigenvalue,damping factor,resolution,electrical abrupt interface

2016-11-08;改回日期2017-06-12。

唐榮江(1992—),男,碩士在讀,主要從事電磁地球物理正、反演方法研究。

國(guó)家自然科學(xué)基金項(xiàng)目(41674078)資助。

This research is financially supported by the National Natural Science Foundation of China (Grant No.41674078).

P631

A

1000-1441(2017)06-0898-07

10.3969/j.issn.1000-1441.2017.06.016