李 賡，曹飛翔，李萬開

(河南理工大學 物理與電子信息學院，焦作 454150)

0 引言

大地電磁測深法是一種天然源頻率域電磁法，其以自然產(chǎn)生的平面電磁波為場源，通過觀測地表正交的電磁場分量獲取地下的電阻率結(jié)構(gòu)[1]。然而，在地形較為復雜的地區(qū)，地面大地電磁法施工效率較低，大面積勘探成本急劇上升。目前，航空大地電磁法具有機動靈活、成本低、效率高的優(yōu)點，能夠在復雜環(huán)境中進行大規(guī)模快速勘探，在國外已被應(yīng)用于油氣勘探、礦產(chǎn)資源勘探等領(lǐng)域[2-3]。

航空大地電磁法通過搭載在直升機上的探測設(shè)備采集垂直磁場信號，同時在地面基站采集水平磁場信號，獲得傾子數(shù)據(jù)。經(jīng)反演[4-5]計算，得到地下介質(zhì)的電性結(jié)構(gòu)。雖然傾子數(shù)據(jù)對電阻率的橫向變化敏感，可以展示橫向電阻率變化引起的異常響應(yīng)，但缺乏背景電阻率的信息。結(jié)合地面大地電磁數(shù)據(jù)和傾子數(shù)據(jù)的反演[6]，可克服只擬合傾子數(shù)據(jù)的缺陷。

近年來，隨著計算機硬件計算能力的提高，有限差分法和有限元法等數(shù)值模擬方法逐步得到應(yīng)用。Sasaki[7]采用有限差分法實現(xiàn)了針對ZTEM傾子響應(yīng)的正演模擬。有限差分法[8-10]在模擬起伏地形或復雜異常體所構(gòu)成的不規(guī)則界面時較困難[11]，而采用有限元方法[12-14]進行起伏地形或不規(guī)則界面模擬時更加容易。

這里實現(xiàn)了二維航空大地電磁有限元正演，分析了起伏地形對航空大地電磁傾子數(shù)據(jù)的影響。為了解決上述傾子數(shù)據(jù)缺乏背景電阻率信息的問題，實現(xiàn)了空中傾子數(shù)據(jù)與地面大地電磁數(shù)據(jù)的反演算法。通過在傾子數(shù)據(jù)的反演目標函數(shù)中加入地面大地電磁數(shù)據(jù)約束項，構(gòu)建了新的反演目標函數(shù)，同時引入加權(quán)因子對地面大地電磁數(shù)據(jù)“約束”項進行加權(quán)。利用非線性共軛梯度(NLCG)反演算法，開展了帶地形的低阻異常體地電模型的傾子反演和傾子與地面大地電磁數(shù)據(jù)反演研究。

1 二維航空大地電磁正演理論

1.1 邊值問題

在航空大地電磁法中，空氣中電阻率一般設(shè)置為108，位移電流的大小跟電導率和頻率成正比。因此，可忽略位移電流對場分布的影響。取時諧因子為e-iωt，于是諧變場的麥克斯韋方程組可表示為：

(1)

由麥克斯韋方程組可進一步推導出

TE極化模式：

(2)

TM極化模式：

(3)

從TE和TM模式極化公式中可以看出，磁場分量只存在于TE模式之中，TE極化模式下二維航空電磁法的控制方程表示為：

(4)

公式(4)可進一步表示成：

(5)

式中：u=Ex;τ=1/iμω;λ=σ。

為了求解式(5)，還須給出邊界條件?，F(xiàn)假設(shè)求解區(qū)域如圖1所示，筆者所研究對象為天然電磁場，場源在距離地面足夠遠的高空中，電磁場以平面波的形式垂直入射地下介質(zhì)。

圖1 TE模式求解區(qū)域Fig.1 Solution domain for TE polarization

TE模式下的邊界條件如下：

1)上邊界AB離地面足夠遠，使異常場在AB上為零，該處的u為“1”。

3)左右邊界AD、BC離局部不均勻體足夠遠，邊界條件為?u/?n。

綜上所述，邊值問題可總結(jié)為式(6)。

u=1 ∈AB

(6)

1.2 有限元分析

與上述邊值問題等價的變分問題為式(7)。

δF(u)=0

(7)

式中：Ω為待求解區(qū)域面積;Γ為Ω的邊界。

有限單元法在求解上述變分問題時，首先將待求解區(qū)域離散為一系列互不重合的四邊形單元。接著選用等參單元并構(gòu)建雙二次插值函數(shù)進行單元積分。圖2展示了母單元(等參單元)和子單元以及其各節(jié)點編號，計算時需要將子單元節(jié)點坐標映射到定義域為[-1,1]2的等參單元上，并用(ξ,η)表示等參變量。設(shè)子單元中心點坐標為(x0,y0)，單元高為b，寬為a，則母單元與子單元之間的坐標映射滿足如下關(guān)系：

圖2 四邊形單元Fig.2 Quadrilateral unit(a)母單元;(b)子單元

(8)

坐標之間微分關(guān)系為式(9)。

(9)

則四邊形母單元的雙二次插值函數(shù)形式如下：

(10)

然后在各離散單元內(nèi)用插值函數(shù)與單元中節(jié)點上未知電場值的乘積作為近似解，將近似解帶入邊值問題后，在每個單元上可得一個復系數(shù)線性方程組。向該復系數(shù)線性方程組中加入求解邊界條件(定解條件)，假設(shè)為u=b。最終可得到一個大型稀疏對稱復系數(shù)線性方程組，該大型稀疏對稱復稀疏線性組表示為式(11)。

(11)

式中:Ke和ue分別為單元系數(shù)矩陣和向量；Ne為單元總個數(shù)，采用直接求解法求解該方程組即可得到求解區(qū)域內(nèi)各節(jié)點電場值組成的向量。

1.3 傾子響應(yīng)計算公式

航空電磁法中傾子響應(yīng)的表示為式(12)。

Hz(r)=Tzx(r,r0)Hx(r0)+

Tzy(r,r0)Hy(r0)

(12)

式中：r為空中垂直磁場分量數(shù)據(jù)采集點位置；r0為地面水平磁場分量數(shù)據(jù)采集點位置；Tzx、Tzy分別為傾子數(shù)據(jù)的兩個分量。

由于磁場分量Hz只存在于TE模式之中，因此，二維航空電磁法中傾子數(shù)據(jù)只存在Tzy分量，其計算公式為式(13)。

(13)

2 反演理論

2.1 反演目標函數(shù)的建立

為實現(xiàn)傾子數(shù)據(jù)與地面大地電磁數(shù)據(jù)的反演，我們在構(gòu)建目標函數(shù)時加入對地面大地電磁數(shù)據(jù)的約束項，如式(14)所示。

(14)

最終該反演目標函數(shù)的完整形式如下：

λmTLTLm

e=d-F(m)

eMT=dMT-FMT(m)

(15)

2.2 反演目標函數(shù)梯度的求取

基于Rodi W[15]的文獻，采用非線性共軛梯度(Nonlinear conjugate gradients,NLCG)算法求解新反演目標函數(shù)的極小值。這里的新目標函數(shù)中，加入了對地面大地電磁數(shù)據(jù)的約束項，此處與經(jīng)典目標函數(shù)不同，如式(15)所示。另外，視電阻率與傾子的計算公式不同，因此在求解正演算子的雅可比矩陣的過程有所不同。

傾子數(shù)據(jù)及正演算子可表示為：

(16)

式中：向量a、b可實現(xiàn)電場值向磁場值的映射；u表示電場分量Ex；雅可比矩陣A是正演算子的偏導數(shù)，即?F(m)/?m，可表示為：

(17)

其中：i=1,2,…,N；j==1,2,…,M。

視電阻率以及正演算子FMT(·)可表示為：

ρa=-Ex/(1/iωμ)(?Ex)/?z)

FMT(m)=lnρa=

(18)

式中：i表示虛數(shù)單位；向量h、k可實現(xiàn)電場值向電場值、磁場值的映射。而該正演算子的偏導數(shù)為雅可比矩陣C,可表示為式(19)。

(19)

其中:i=1,2…,NMT;j=1,2…,M。

3 算例

3.1 數(shù)值模擬算例

3.1.1 正確性驗證

為了驗證本文有限元算法的正確性，設(shè)計了如3所示的二維低阻塊體模型。該模型的背景電阻率為100 Ω·m，異常體的電阻率為10 Ω·m，其頂層距離地表1 000 m，橫向?qū)挾葹? 000 m，垂向高度為 2 000 m。空中測線位于地表以上100 m處(圖3中紅色虛線)，頻率選擇為0.01 Hz、1 Hz和100 Hz，分別計算了傾子實部和虛部，并將計算結(jié)果與MARE2DEM程序[16]計算結(jié)果進行了對比如圖4、圖5所示。

圖3 低阻塊體模型Fig.3 Low resistivity block model

圖4 傾子實部對比曲線圖Fig.4 The real part of tipper comparison curve

圖5 傾子虛部對比曲線圖Fig.5 The imaginary part of tipper comparison curve

圖6 傾子實部相對誤差曲線圖Fig.6 The real part of tipper relative error curve

圖7 傾子虛部相對誤差曲線圖Fig.7 The imaginary part of tipper relative error curve

圖8 起伏地形模型示意圖Fig.8 The schematic diagram of undulating terrain model

圖9 傾子響應(yīng)實部Fig.9 The real part of tipper

圖6、圖7分別展示了各不同頻率傾子實虛部相對誤差曲線，從圖6、圖7中可以看出各分量相對誤差均小于0.7%，說明本文數(shù)值模擬算法的正確性。

3.1.2 起伏地形

圖10 傾子響應(yīng)虛部Fig.10 The imaginary part of tipper

圖11 異常體反演模型Fig.11 Anomaly model for inversion

為了分析傾子數(shù)據(jù)在起伏地形環(huán)境下的響應(yīng)特征，將圖8所示的起伏地形正演模型(背景電阻率值為500 Ω·m的均勻半空間)中的山谷和山脊地形，設(shè)計為地形范圍位于-3 000 m≤y≤-1 000 m和1 000 m≤y≤3 000 m內(nèi)，最大地形為500 m?？罩袦y線位于地形以上100 m，如圖8中紅色虛線所示。在-4 000 m≤y≤4 000 m范圍內(nèi)每間隔50 m采集磁場垂直分量，在地表y=-5 000 m處采集磁場水平分量。選擇30 Hz、103.923 Hz和360 Hz頻率計算傾子響應(yīng)值。圖9和圖10分別展示了起伏地形各頻率傾子響應(yīng)值的實虛部。

通過分析起伏地形傾子響應(yīng)可知，起伏地形的存在使傾子實虛部響應(yīng)值曲線偏離于零值(平地形傾子響應(yīng)為零)，這種地形效應(yīng)主要出現(xiàn)在山谷和山脊地形處，且隨著頻率的增加，傾子響應(yīng)曲線的偏離程度逐漸增大。這是由于高頻段信號攜帶的淺部介質(zhì)的電性結(jié)構(gòu)信息較多，當?shù)乇硖庪娦越Y(jié)構(gòu)發(fā)生變化時，高頻段信號也會隨之發(fā)生明顯變化，而低頻段信號反映的是深部介質(zhì)的電性結(jié)構(gòu)，地表處的電性結(jié)構(gòu)變化對其影響較小。

3.2 反演算例

我們對圖11中帶異常體的起伏地形模型進行反演。模型的背景電阻率為500 Ω·m,在山脊和山谷下方分別嵌入一個電阻率為10 Ω·m的異常體，其高度為500 m，寬度均為800 m，山谷地形下的異常體頂面距離其地形最低處750 m，山脊地形下的異常體頂面距離其地形最高處1 250 m。圖11中，航空測線位于地形以上100 m。

圖12 500 Ω·m均勻半空間初始模型航空大地電磁反演結(jié)果Fig.12 Inversion result of Airborne Magnetotelluric with 500 Ω·m uniform half-space

圖13 2 000 Ω·m均勻半空間初始模型航空大地電磁反演結(jié)果Fig.13 Inversion result of Airborne Magnetotelluric with 2 000 Ω·m uniform half-space

在-4 000 m≤y≤4 000 m范圍內(nèi)按照20 m間隔共選取401個航空測點，采集空中磁場垂直分量。在地表y=-5 000 m處觀測磁場水平分量。在30 Hz～360 Hz頻率范圍內(nèi)對數(shù)等間隔選擇11個計算頻率，共生成401×11×2(傾子實部和虛部)=8 822個傾子數(shù)據(jù)。同時為了進行地面大地電磁數(shù)據(jù)的聯(lián)合反演，在-4 000 m≤y≤4 000 m范圍內(nèi)，共放置了12個地面大地電磁數(shù)據(jù)采集站。在1 Hz～1 000 Hz頻率范圍內(nèi)對數(shù)等間隔選擇13個計算頻率值，共生成12×13×2(視電阻率和阻抗相位)=312個地面大地電磁數(shù)據(jù)。分別向數(shù)值模擬計算得到的傾子數(shù)據(jù)和地面大地電磁數(shù)據(jù)中加入3%的隨機噪聲充當實際觀測數(shù)據(jù)。

3.2.1 帶地形的傾子反演

均勻半空間模型的傾子為零，傾子的梯度不為零。在進行帶地形的傾子反演時，分別將初始模型設(shè)置為電阻率值為500 Ω·m和2 000 Ω·m的均勻半空間，其對應(yīng)的傾子數(shù)據(jù)反演結(jié)果分別如圖12和圖13所示，由圖12可知，當反演初始模型電阻率值與真實背景電阻率值相同時，由于地形的影響，傾子數(shù)據(jù)反演結(jié)果中異常體的位置相對于真實位置整體偏上，異常體的電阻率值也比真實值偏大。

圖14 500 Ω·m均勻半空間初始模型聯(lián)合反演反演結(jié)果Fig.14 Inversion result of Joint inversion with 500 Ω·m uniform half-space

圖15 2 000 Ω·m均勻半空間初始模型聯(lián)合反演反演結(jié)果Fig.15 Inversion result of Joint inversion with 2 000 Ω·m uniform half-space

由圖13可知，當反演初始模型電阻率值與真實背景電阻率值不同時，傾子數(shù)據(jù)的反演結(jié)果與真實電阻率模型相差較大，反演結(jié)果中異常體位置、大小未得到有效反映。這是由于傾子數(shù)據(jù)雖然可以有效表示電阻率的相對變化，但是對電阻率的絕對值約束力不足，當初始模型與真實模型相差較大時，單獨的傾子數(shù)據(jù)反演無法得到一個可靠的反演結(jié)果。

3.2.2 傾子與地面大地電磁數(shù)據(jù)反演

為驗證本文反演算法的必要性，我們針對一組加入了相同噪聲的傾子和地面大地電磁合成數(shù)據(jù)，采用不同的反演初始模型分別進行聯(lián)合反演。在進行傾子與地面大地電磁數(shù)據(jù)反演時，分別將反演初始模型選擇為電阻率為500 Ω·m和2 000 Ω·m的均勻半空間，其反演結(jié)果分別如圖14和圖15所示。

當反演初始模型電阻率值和真實背景電阻率相同時，向傾子數(shù)據(jù)中加入少量地面大地電磁數(shù)據(jù)進行反演后，異常體的橫向邊界可以被清晰的反映出來，同時異常體的反演電阻率值得到了改善，反演效果得到了一定程度的提高。

圖13是初始反演背景電阻率值為2 000 Ω·m的傾子反演結(jié)果，分析圖13和圖15可知，加入少量的地面大地電磁數(shù)據(jù)進行反演，反演結(jié)果得到了明顯改善。不但兩個低阻異常體的邊界得到了有效刻畫，而且其反演電阻率值也與異常體真實電阻率值較接近。同時對比圖14與圖15可以發(fā)現(xiàn)，當初始模型電阻率與真實模型背景電阻率值差異較大時，少量地面大地電磁數(shù)據(jù)與傾子數(shù)據(jù)的反演結(jié)果并沒有發(fā)生巨大改變，這說明通過加入少量的地面大地電磁數(shù)據(jù)進行反演，可有效地降低傾子數(shù)據(jù)反演結(jié)果對初始模型的依賴。

4 結(jié)論

首先采用四邊形網(wǎng)格模擬起伏地形，使用有限元方法對起伏地形均勻半空間模型進行數(shù)值模擬，分析和總結(jié)了起伏地形對航空大地電磁傾子數(shù)據(jù)的影響。接著對起伏地形異常體模型進行了傾子反演和傾子與少量地面大地電磁數(shù)據(jù)反演。最終結(jié)論如下：

1)起伏地形的存在使得傾子實虛部響應(yīng)值曲線不同程度的偏離于零值，且隨著計算頻率的增加，傾子響應(yīng)曲線的偏離程度逐漸增大。這種地形效應(yīng)嚴重影響了航空大地電磁傾子響應(yīng)曲線的形態(tài)，起伏地形不可忽視，會直接影響到航空大地電磁數(shù)據(jù)資料的處理解釋結(jié)果。

2)利用傾子數(shù)據(jù)進行反演時，由于傾子數(shù)據(jù)對電阻率的絕對值約束力不足，反演結(jié)果受初始模型影響嚴重，當初始模型與真實模型相差很大時，單獨的傾子數(shù)據(jù)反演無法得到可靠的反演結(jié)果；通過加入少量的地面大地電磁數(shù)據(jù)進行反演，可增加觀測完備性，從而降低多解性，最終可有效降低傾子數(shù)據(jù)反演結(jié)果對初始模型的依賴。