張雙全，李桐林

吉林大學 地球探測科學與技術(shù)學院，長春 130026

0 引言

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

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

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

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

1 方法理論

1.1 正演算法

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

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

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

TE模式：

(1)

TM模式：

(2)

從以上兩個公式可見，任意各向異性情況下，

TE模式：

(3)

TM模式：

(4)

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

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

1.2 反演算法

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

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

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

(5)

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

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

(6)

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

[dobs-dk+Jmk-Jmref]

(7)

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

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

模型1

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

計算時將模型剖分為72×24的剖分單元(垂向24個剖分單元，不用包含空氣層)，測點數(shù)為72個(即地表的網(wǎng)格節(jié)點處)，采用25個周期點(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é)果來看(圖2)，y軸方向上成功反演出了異常體的水平位置和深度。在實際異常區(qū)域內(nèi)值大約為160 Ω·m，略小于真實值200 Ω·m，這是因為MT平面波的性質(zhì)，即極化電場和磁場在水平方向衰減，且在高阻中衰減緩慢。z軸方向上反演的ρzz與真實模型不一致，而是與背景空間接近。

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

模型2

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

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

模型3

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

計算時將模型剖分為52×24的剖分單元(垂向24個剖分單元，不用包含空氣層)，測點數(shù)為52個(即地表的網(wǎng)格節(jié)點處)，采用25個周期點(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é)果來看(圖5)，y軸方向上成功反演出了兩個異常體的水平位置和深度。ρyy在實際異常區(qū)域內(nèi)值大約為150 Ω·m，略小于真實值200 Ω·m。z軸方向上反演的ρzz與真實模型不一致，而是與背景空間接近。

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

模型4

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

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

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

3 結(jié)論

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

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