二維大地電磁主軸各向異性理論模型反演

張雙全，李桐林

吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

大地電磁法(maganetotelluric method，MT)利用天然存在的大地電磁場(chǎng)，能反映出地殼和上地幔范圍內(nèi)的地下電性結(jié)構(gòu)，具有其他勘探方法無(wú)法比擬的優(yōu)勢(shì)。電性各向異性是指電導(dǎo)率隨方位變化的現(xiàn)象。Wannamaker[1]詳細(xì)論述了電性各向異性在地電模型和地球動(dòng)力學(xué)解釋中的重要作用，還說(shuō)明了各向異性介質(zhì)普遍存在于地球內(nèi)部。但在實(shí)際應(yīng)用的大地電磁數(shù)據(jù)采集和分析解釋流程中，通常直接簡(jiǎn)單地將地下介質(zhì)設(shè)定為各向同性情況而忽略電性各向異性現(xiàn)象，這樣很容易漏掉甚至得到不真實(shí)的地質(zhì)信息[2]。

20世紀(jì)60年代末對(duì)電導(dǎo)率各向異性的理論研究逐漸發(fā)展起來(lái)。Reddy et al.[3]首次對(duì)二維電性各向異性偏微分方程用有限元法做了數(shù)值計(jì)算，對(duì)后世的理論研究留下了尤為深遠(yuǎn)的影響。進(jìn)入20世紀(jì)90年代在計(jì)算機(jī)技術(shù)的輔助下，電磁數(shù)值模擬技術(shù)飛速發(fā)展，電各向異性研究也取得重大突破[4]，有限元法和有限差分法被廣泛用于二維大地電磁各向異性的正演模擬中[5]。其中Pek J et al.[6]

在Reddy et al.[3]基礎(chǔ)上，將前人Haak[7]，Brewitt--

中國(guó)關(guān)于電性各向異性的研究雖然晚于國(guó)外，但也取得了許多成果。Li[11]公布了二維各向異性電磁法有限元算法，并在后來(lái)與Pek J用自適應(yīng)的網(wǎng)格剖分法將算法做了改進(jìn)[12]，是二維各向異性有限元電磁模擬中影響較為重大的研究成果。楊長(zhǎng)福等[13]在Pek J的理論基礎(chǔ)上，將電導(dǎo)率簡(jiǎn)化為對(duì)稱各向異性介質(zhì)，反演了甘肅天祝永登一帶的大地電磁資料數(shù)據(jù)并取得了不錯(cuò)的效果。胡祥云等[14]將Pek J的研究方法應(yīng)用到新疆某地的大地電磁資料處理中，做了二維正演擬合解釋，驗(yàn)證了算法的實(shí)用性和理論的正確性。由于各向異性反演較為復(fù)雜，單純利用大地電磁各向異性的應(yīng)用案例較少，筆者編寫(xiě)了Matlab計(jì)算程序，并反演計(jì)算了多個(gè)二維模型，為更廣泛地應(yīng)用大地電磁各向異性提供了理論依據(jù)。本研究存在一定局限性，只是考慮了主軸各向異性的情況，忽視了各主軸與構(gòu)造的夾角，而且反演的是較為簡(jiǎn)單的理論模型，沒(méi)有研究大地電磁實(shí)測(cè)數(shù)據(jù)資料。介于目前國(guó)際上還未出現(xiàn)成熟的、可同時(shí)反演3個(gè)主軸電阻率和旋轉(zhuǎn)角的算法，且各向異性識(shí)別方面的研究也不夠完善，導(dǎo)致電性各向異性在實(shí)際生產(chǎn)資料處理中應(yīng)用不多，因此不考慮旋轉(zhuǎn)角，只計(jì)算主軸電阻率，不失為一個(gè)較好的突破口，也有一定的應(yīng)用價(jià)值，能為今后處理解釋復(fù)雜的大地電磁資料的電性各向異性現(xiàn)象提供思路和方法。

1 方法理論

1.1 正演算法

本文使用的二維大地電磁各向異性有限差分正演理論來(lái)自Pek et al.[6]的研究方法，已有許多學(xué)者驗(yàn)證了其正確性。

對(duì)任意各向異性異常體，假設(shè)在笛卡爾坐標(biāo)系中，其走向平行于x軸，y軸垂直于x軸，z軸垂直于xy平面且向下為正向。地表水平且對(duì)應(yīng)z=0，地表上空是空氣層，來(lái)自高空的平面波垂直地表向下傳播，有：

它可分解為3個(gè)主軸方向的電導(dǎo)率和一個(gè)由3個(gè)Euler’s空間旋轉(zhuǎn)角αL、αD和αS組成的旋轉(zhuǎn)矩陣。將其帶入麥克斯韋方程組，忽略位移電流后，可得到以下一組關(guān)于電場(chǎng)分量Ex和磁場(chǎng)水平分量Hx的二階偏微分方程(嚴(yán)格來(lái)說(shuō)在各向異性情況下，電場(chǎng)和磁場(chǎng)分量耦合，已不存在“純粹”的TE模式與TM模式，不過(guò)為表述方便，下文依然將電場(chǎng)沿走向的E極化模式稱作“TE模式”，將磁場(chǎng)沿走向的H極化模式稱作“TM模式”)：

TE模式：

(1)

TM模式：

(2)

從以上兩個(gè)公式可見(jiàn)，任意各向異性情況下，

TE模式：

(3)

TM模式：

(4)

從公式(3)和(4)可以看出，電場(chǎng)分量Ex和磁場(chǎng)分量Hx已成功解耦，為單獨(dú)反演電場(chǎng)分量和磁場(chǎng)分量提供了提論依據(jù)。TE模式只反映σxx且在形式上公式(3)與二維各向同性介質(zhì)相同，只是在物理意義上公式(3)中反映的是沿走向方向的電導(dǎo)率值，這意味著此時(shí)二維各向同性介質(zhì)的數(shù)值模擬方法甚至程序可直接利用，還可推廣至反演[13]。二維大地電磁各向同性介質(zhì)的正反演已非常成熟，因此本文不再研究二維各向異性TE模式的計(jì)算。相比各向同性介質(zhì)，二維各向異性介質(zhì)情況下的TM模式就復(fù)雜一些，公式(4)可以看出TM模式依然同時(shí)與σyy和σzz兩個(gè)電導(dǎo)率值相關(guān)，且σyy≠σzz，各向同性介質(zhì)下的計(jì)算程序無(wú)法直接使用。

相應(yīng)地， TE模式下每個(gè)節(jié)點(diǎn)列方程時(shí)的10個(gè)系數(shù)只剩下和電場(chǎng)分量有關(guān)的5個(gè)，TM模式下的14個(gè)系數(shù)也只剩下和磁場(chǎng)分量有關(guān)的5個(gè)。

1.2 反演算法

針對(duì)反演問(wèn)題已有多種數(shù)值優(yōu)化技術(shù)，如Occam、非線性共軛梯度、高斯牛頓、擬牛頓法和最小二乘法等[4]。本文采用最小二乘法，使用基于先驗(yàn)?zāi)Ｐ偷淖罟饣Ｐ图s束。

反演的總目標(biāo)函數(shù)可以表示為：

Pα(m)=‖[dobs-F(m)]‖2+α‖Wm(m-mref)‖2

(5)

式中：Pα(m)表示總目標(biāo)函數(shù)；α是正則化因子；F表示正演響應(yīng)，Wm是光滑度矩陣，或稱模型權(quán)系數(shù)矩陣；mref是先驗(yàn)?zāi)Ｐ汀?/p>

Pα(mk+1)=‖[dobs-dk-JkΔm‖2+α‖Wm(mk-mref)‖2

(6)

上述函數(shù)對(duì)Δm求導(dǎo)并令它等于0，寫(xiě)成迭代形式為:

[dobs-dk+Jmk-Jmref]

(7)

J是雅可比矩陣，也稱靈敏度矩陣，代表每個(gè)數(shù)據(jù)對(duì)模型參數(shù)的偏微分。雅可比矩陣可用差分或伴隨陣等方法計(jì)算得到。對(duì)公式(7)進(jìn)行反演迭代，將計(jì)算得到的優(yōu)化模型改變量不斷更新，重復(fù)計(jì)算目標(biāo)函數(shù)，直至目標(biāo)函數(shù)減小到符合精度要求即可。

2 二維主軸各向異性模型計(jì)算與分析

模型1

圖1 理論模型1Fig.1 Theoretical model 1

計(jì)算時(shí)將模型剖分為72×24的剖分單元(垂向24個(gè)剖分單元，不用包含空氣層)，測(cè)點(diǎn)數(shù)為72個(gè)(即地表的網(wǎng)格節(jié)點(diǎn)處)，采用25個(gè)周期點(diǎn)(0.01，0.03，0.05，0.07，0.1，0.3，0.5，0.7，1，3，5，7，10，30，50，70，100，300，500，700，1 000，1 300，1 500，1 700，2 000，單位s)。反演初始模型設(shè)置為與介質(zhì)1完全相同的各向同性均勻半空間。

從反演結(jié)果來(lái)看(圖2)，y軸方向上成功反演出了異常體的水平位置和深度。在實(shí)際異常區(qū)域內(nèi)值大約為160 Ω·m，略小于真實(shí)值200 Ω·m，這是因?yàn)镸T平面波的性質(zhì)，即極化電場(chǎng)和磁場(chǎng)在水平方向衰減，且在高阻中衰減緩慢。z軸方向上反演的ρzz與真實(shí)模型不一致，而是與背景空間接近。

(a)y軸方向電阻率反演結(jié)果；(b)z軸方向電阻率反演結(jié)果。圖2 理論模型1的反演結(jié)果Fig.2 Inversion results for theoretical model 1

模型2

從反演結(jié)果來(lái)看(圖3)，與圖2類似，y軸方向上的效果好于z軸。y軸方向上除了較好地反演出了異常體的位置，還成功反演出了異常體的真實(shí)值為20 Ω·m，而z軸方向上則沒(méi)有取得較好的效果。

(a)y軸方向電阻率反演結(jié)果；(b)z軸方向電阻率反演結(jié)果。圖3 理論模型2的反演結(jié)果Fig.3 Inversion results for theoretical model 2

模型3

圖4 理論模型3Fig.4 Theoretical model 3

計(jì)算時(shí)將模型剖分為52×24的剖分單元(垂向24個(gè)剖分單元，不用包含空氣層)，測(cè)點(diǎn)數(shù)為52個(gè)(即地表的網(wǎng)格節(jié)點(diǎn)處)，采用25個(gè)周期點(diǎn)(0.01，0.03，0.05，0.07，0.1，0.3，0.5，0.7，1，3，5，7，10，30，50，70，100，300，500，700，1 000，1 300，1 500，1 700，2 000，單位s)。反演初始模型設(shè)置為與介質(zhì)1完全相同的各向同性均勻半空間。

從反演結(jié)果來(lái)看(圖5)，y軸方向上成功反演出了兩個(gè)異常體的水平位置和深度。ρyy在實(shí)際異常區(qū)域內(nèi)值大約為150 Ω·m，略小于真實(shí)值200 Ω·m。z軸方向上反演的ρzz與真實(shí)模型不一致，而是與背景空間接近。

(a)y軸方向電阻率反演結(jié)果；(b)z軸方向電阻率反演結(jié)果。圖5 理論模型3的反演結(jié)果Fig.5 Inversion results for theoretical model 3

模型4

從反演結(jié)果來(lái)看(圖6)，與模型3類似，y軸方向上成功反演出了兩個(gè)異常體的水平位置和深度；z軸方向上反演的ρzz與真實(shí)模型不一致，而是與背景空間接近。右側(cè)高阻模型的ρyy在實(shí)際異常區(qū)域內(nèi)的值約為140 Ω·m，小于真實(shí)值200 Ω·m，且小于模型3中的150 Ω·m，說(shuō)明右側(cè)的高阻體受到了左側(cè)低阻體的影響；左側(cè)低阻模型的ρyy在實(shí)際異常區(qū)域內(nèi)的值約為20 Ω·m，與真實(shí)值相等。

(a)y軸方向電阻率反演結(jié)果；(b)z軸方向電阻率反演結(jié)果。圖6 理論模型4的反演結(jié)果Fig.6 Inversion results for theoretical model 4

綜合以上4個(gè)反演結(jié)果，可見(jiàn)MT響應(yīng)對(duì)低阻異常更為敏感。相較高阻，反演能更有效地恢復(fù)低阻電阻率；相較縱向，水平方向異常體的范圍圈定要更準(zhǔn)確些。同時(shí)這4個(gè)模型還說(shuō)明了在二維各向異性反演中恢復(fù)ρzz幾乎是不可能實(shí)現(xiàn)的，這是由于其非固有唯一性導(dǎo)致的。Yin[15]證明了主軸各向異性情況下ρzz不可能反演恢復(fù)；Chen[16]在二維各向異性反演中得到了相同的結(jié)論。

3 結(jié)論

本文編程實(shí)現(xiàn)了二維大地電磁各向異性“TM”模式下理論模型的最小二乘反演。反演計(jì)算了單個(gè)異常體模型和多個(gè)異常體模型，在y軸方向上均得到了較好的反演結(jié)果：位置圈定上，異常體的范圍與理論模型非常接近，水平方向上的范圍極為準(zhǔn)確，縱向上稍差一些；電阻率上，高阻體會(huì)受到低阻體的影響，低阻的電阻率恢復(fù)相較高阻更加準(zhǔn)確。反演結(jié)果同時(shí)驗(yàn)證了前人關(guān)于z軸方向上的電阻率由于ρzz不敏感，無(wú)法成功恢復(fù)的研究結(jié)論。

致謝感謝捷克科學(xué)院地球物理研究所Josef Pek教授對(duì)正演算法的提供。