李志強(qiáng)，孫洋，譚捍東，張承客

(1.江西省交通科學(xué)研究院，江西 南昌 330200； 2.中國地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院，北京 100083)

0 引言

Z軸傾子電磁法(Z-axis tipper electromagnetic，ZTEM)是將天然場勘探深度大和航空測量工作效率高的特點(diǎn)有效地結(jié)合起來，基于音頻磁場法AFMAG[1-2]發(fā)展而來的一種新型頻率域航空電磁法[3]。該方法利用航空測量與地面基站觀測的野外勘探方式，接收天然場源下25～720 Hz頻段的磁場數(shù)據(jù)，并通過傾子T[4]建立磁場之間的聯(lián)系，可以快速獲得起伏地形下三維電導(dǎo)率構(gòu)造信息。相比人工源航空電磁法，ZTEM具有勘探深度大、經(jīng)濟(jì)、高效等優(yōu)點(diǎn)。

由于ZTEM是近年來才提出的一種航空電磁法，目前國外一些學(xué)者在從事這種方法的二維與三維正反演研究，并且應(yīng)用到實(shí)例當(dāng)中檢驗(yàn)算法的有效性，目的是想利用頻率域航空電磁法探測深部的地質(zhì)構(gòu)造和礦體，并且已取得較為豐富的成果[5-8]，實(shí)踐應(yīng)用表明ZTEM方法是實(shí)現(xiàn)經(jīng)濟(jì)快速探礦的理想手段。國內(nèi)方面，王濤等[9]為探索該方法的參數(shù)響應(yīng)特征，采用有限差分算法進(jìn)行了ZTEM三維正演研究，為ZTEM的反演研究奠定了較好的基礎(chǔ)。趙叢等[10]對比介紹了ZTEM方法和MT方法的特點(diǎn)，并提出利用這兩種天然場源電磁法進(jìn)行點(diǎn)面結(jié)合的聯(lián)合深部礦產(chǎn)資源勘探思路。李志強(qiáng)等[11]采用三維數(shù)據(jù)空間OCCAM算法進(jìn)行了ZTEM傾子資料的合成算例反演研究，并且與MT阻抗反演結(jié)果對比檢驗(yàn)了算法的有效性，同時(shí)表明ZTEM參數(shù)反演的橫向約束能力突出。而針對ZTEM空心感應(yīng)磁傳感器的感應(yīng)線圈具有源阻抗大的特點(diǎn)，王言章等[12]對ZTEM磁傳感器調(diào)理電路做了低噪聲優(yōu)化設(shè)計(jì)。此外，許智博等[13]基于有限差分正演還實(shí)現(xiàn)了二維ZTEM非線性共軛梯度(NLCG)反演算法。

考慮地形因素的影響，國外學(xué)者Sasaki等[14]證明了ZTEM三維反演時(shí)地形效應(yīng)在較高頻時(shí)更顯著的特點(diǎn)，Legault等[15]和Sattel等[16]分別對比分析了起伏地形下二維、三維ZTEM和MT單獨(dú)與聯(lián)合反演的效果，說明了聯(lián)合反演中可利用MT反演結(jié)果來校正ZTEM反演時(shí)的背景電阻率，而三維聯(lián)合反演更能真實(shí)恢復(fù)出地下電導(dǎo)率結(jié)構(gòu)。當(dāng)前國內(nèi)研究地形對航空電磁的影響主要集中在主動(dòng)源的正反演數(shù)值模擬中，張博等[17]應(yīng)用非結(jié)構(gòu)網(wǎng)格矢量有限元法模擬分析了起伏地表下的三維頻域/時(shí)域航空電磁響應(yīng)特征，為航空電磁地形效應(yīng)的有效識別提供了可能；王衛(wèi)平等[18]采用有限差分法分析了頻率域航空電磁法的地形影響，并采用地形校正方法對典型地形地電模型進(jìn)行了校正，提高了異常響應(yīng)特征的解釋效果；李文奔[19]開發(fā)實(shí)現(xiàn)了起伏地表下的頻率域二維航空電磁高斯-牛頓反演算法，通過模型算例表明忽略地形反演會(huì)造成虛假異常的產(chǎn)生。國內(nèi)學(xué)者為了適應(yīng)地形起伏的實(shí)際情況，也對天然場源的大地電磁法(MT)進(jìn)行了帶地形的二維、三維正反演研究[20-21]。鑒于航空電磁法在起伏地形區(qū)域的勘探優(yōu)勢明顯，而實(shí)際野外探測作業(yè)中地形的變化影響又不容忽視，本文基于ZTEM三維數(shù)據(jù)空間OCCAM反演算法[11]，主要研究了帶地形的ZTEM傾子數(shù)據(jù)三維正反演算法，通過合成算例計(jì)算和分析純地形的ZTEM三維響應(yīng)(傾子)特征，對典型地電模型ZTEM帶地形與不帶地形的反演算例進(jìn)行了探討分析，檢驗(yàn)了帶地形的ZTEM數(shù)據(jù)空間OCCAM反演算法的有效性與可靠性。

1 起伏地形下ZTEM三維正演

1.1 傾子資料整理

ZTEM方法的傾子數(shù)據(jù)是將地面以上空氣層中不同位置處r的航空測量垂直磁場與地表某一固定參考基站處r0的兩個(gè)水平磁場聯(lián)系起來[11]，其關(guān)系表示如下：

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

(1)

式中：Tzx和Tzy表示傾子分量。

(2)

解方程(2)即可得到傾子T的關(guān)系式：

(3)

1.2 傾子響應(yīng)的計(jì)算

在音頻段的ZTEM勘探方法中，位移電流可忽略，假定地下介質(zhì)磁導(dǎo)率近似為自由空間中的磁導(dǎo)率μ0，時(shí)諧因子取為e-iωt，則ZTEM方法電強(qiáng)度E和磁場強(qiáng)度H所滿足的Maxwell方程組的微分關(guān)系式可表示為

(4)

(5)

式中：ω是角頻率，μ0是磁導(dǎo)率(真空中)，σ是電導(dǎo)率。

KE=S,

(6)

式中：K為矩陣元素與網(wǎng)格單元電阻率及尺寸有關(guān)的對稱大型稀疏系數(shù)矩陣；E表示需要求解的電場三分量列向量；S為列向量，與場源及邊界值有關(guān)。在場源SX極化模式下，方程左端的電場值記為E(1)，右端列向量記為S(1)；在SY極化模式下，方程左端的電場值記為E(2)，右端列向量記為S(2)，則式(6)可分解為兩個(gè)正演方程：

KE(1)=S(1),KE(2)=S(2)。

(7)

圖1 編號為(i，j，k)的網(wǎng)格單元Fig.1 The (i，j，k) grid cell

研究區(qū)域的電場值E(1)和E(2)就通過QMR法求解這兩個(gè)正演方程得到。

(8)

式中：hG表示區(qū)分有兩個(gè)對應(yīng)兩種極化方式下磁場與電場轉(zhuǎn)化關(guān)系式G向量,h=1或2。正演模擬時(shí)，對于起伏地表水平磁場觀測基站所處單元為(i0,j0,ki0,j0)上表面中心處Hj的兩個(gè)分量Hjx、Hjy,對應(yīng)方程(8)的關(guān)系式[25-26]分別為：

(9a)

(9b)

空氣層中某網(wǎng)格單元在(i,j,kair)面中心上的垂直磁場Hjz可表示為：

(Ey(i+1,j,kair)-Ey(i,j,kair))/Δx(i)]。

(9c)

1.3 正演精度的驗(yàn)證

選取二維四棱臺狀山峰純地形模型來試算，以此驗(yàn)證帶地形的正演計(jì)算結(jié)果的準(zhǔn)確性。如圖2所示，走向?yàn)閤方向，四棱臺狀山峰在yz方向上頂面距地表高160 m，上頂邊長800 m，下底邊長3 200 m，地下介質(zhì)與空氣層電阻率分別設(shè)為100 Ω·m和108Ω·m。選取(10 100 m，0 m)地表固定點(diǎn)作為兩個(gè)水平磁場分量取值位置，垂直磁場的取值高度設(shè)置在離四棱臺狀山峰頂面高100 m的空氣層中。

用上述帶地形的三維有限差分正演算法計(jì)算出純地形模型的數(shù)值解，將其結(jié)果與二維有限差分法計(jì)算得到的傾子響應(yīng)進(jìn)行對比，圖3a、b表示兩種數(shù)值模擬算法當(dāng)頻率在25 Hz與200 Hz時(shí)ZTEM異常響應(yīng)Tzy分量的實(shí)、虛部對比結(jié)果，其中“▽”和“○”分別表示頻率在25 Hz和200 Hz時(shí)的三維數(shù)值解，“*”和“×”分別表示頻率在25 Hz和200 Hz時(shí)的二維數(shù)值解。圖3中二維與三維數(shù)值解對比反映出傾子響應(yīng)的實(shí)、虛部都能較好的吻合，證明了所開發(fā)的ZTEM帶地形的三維正演算法準(zhǔn)確可靠。

圖2 ZTEM正演模型Fig.2 The model of ZTEM Forward modeling

2 三維反演算法

本文帶地形的反演算法是基于ZTEM三維數(shù)據(jù)空間OCCAM反演代碼的基礎(chǔ)上所實(shí)現(xiàn)的，實(shí)際反演問題中，假定觀測數(shù)據(jù)總個(gè)數(shù)為N，模型單元總個(gè)數(shù)為M，將目標(biāo)函數(shù)[26]定義為：

(10)

式中:m為待更新的電阻率模型，m0為帶地形的初始模型，Cm和Cd分別是模型協(xié)方差矩陣和數(shù)據(jù)協(xié)方差矩陣，λ-1表示拉格朗日乘子，d表示觀測數(shù)據(jù)，F(xiàn)[m]是帶地形的模型響應(yīng)，X*為期望擬合差水平。

對于第k+1次迭代，將模型空間OCCAM的迭代公式[22]作以下數(shù)學(xué)變換：

(11)

(12)

圖3 二維山峰純地形模型傾子響應(yīng)的二維、三維計(jì)算結(jié)果對比Fig.3 The tipper response comparison between the 2D finite difference modeling results and the 3D finite difference modeling results generated from 2D peak terrain model

帶地形的ZTEM數(shù)據(jù)空間OCCAM反演問題中，為實(shí)現(xiàn)模型的更新迭代，需要通過計(jì)算得到的雅克比矩陣Jk獲得每次模型更新所需的模型修正量，而最為關(guān)鍵的是計(jì)算正演響應(yīng)F和求解對應(yīng)的Jk。傾子用于反演的第k個(gè)模型參數(shù)的偏導(dǎo)數(shù)Jk=(?F/?m)k可相應(yīng)表示為Jk=(?T/?m)k[10]，將式(3)對模型電導(dǎo)率參數(shù)σk求偏導(dǎo)數(shù)，有:

圖4 數(shù)據(jù)空間OCCAM反演流程Fig.4 The flow chart of data-space OCCAM inversion

(13)

由于正演模擬采用的是電場方程，首先將電場正演方程(6)KE=S對模型電導(dǎo)率σk求偏導(dǎo)，考慮S與σk之間的弱相關(guān)性，可以將?S/?σk略去，有

?E/?σk=K-1(-?K/?σkE)。

(14)

(15)

接著將關(guān)系式(15)代入式(13)，則第k個(gè)模型電導(dǎo)率參數(shù)σk的擾動(dòng)對第j個(gè)測點(diǎn)的改變量?Tzxj/?σk與?Tzyj/?σk就寫成如下形式：

式中：

(16)

在ZTEM正演中，差分系數(shù)對稱矩陣K滿足(K-1)T=K-1，故上式中?Tzxj/?σk和?Tzyj/?σk可寫成:

?Tzxj/?σk=(-?K/?σkE1)TK-1(1Gjzx)+

(-?K/?σkE2)TK-1(2Gjzx),

(17)

?Tzyj/?σk=(-?K/?σkE1)TK-1(1Gjzy)+

(-?K/?σkE2)TK-1(2Gjzy)。

(18)

從式(17)、(18)可看出，先把1Gjzx、2Gjzx計(jì)算出后再求解

K·v1=1Gjzx,

(19)

K·v2=2Gjzx,

(20)

可得到向量v1=K-1(1Gjzx)和向量v2=K-1(2Gjzx)，同時(shí)解兩次正演方程獲得電場值E(1)、E(2)后代回式(17)，即可求得某個(gè)測點(diǎn)位置的傾子分量Tzx的偏導(dǎo)數(shù)?Tzxj/?σk。同樣方法可以計(jì)算求出?Tzyj/?σk。

3 合成算例分析

3.1 三維純地形的異常特征

設(shè)計(jì)一個(gè)正四棱臺狀的山峰地形來分析其ZTEM正演傾子平面響應(yīng)特征。如圖5所示，棱臺狀山峰的頂面相對于地面有0.45 km，上頂面邊長有0.5 km，下底邊長有2.5 km，空氣層電阻率為108Ω·m， 地下介質(zhì)的電阻率為100 Ω·m， 選取地表固定點(diǎn)(5 750 m，5 750 m，0 m)作為兩個(gè)水平磁場分量取值位置，垂直磁場的取值高度設(shè)置在離山峰高100 m處。x、y、z方向網(wǎng)格單元剖分?jǐn)?shù)分別為31、31、31(包含空氣層數(shù)為14)。

圖5 山峰地形模型在xy方向的模型視圖(a)、xz方向的模型視圖(b)及xyz方向的模型視圖(c)Fig.5 Schematic diagram of peak terrain model in the view in xy direction(a), in the view in xz direction, and the view in xyzdirection

圖6展示了頻率為50 Hz山峰純地形的ZTEM三維傾子響應(yīng)，圖中白色實(shí)線表示山峰下底邊界。由圖不難發(fā)現(xiàn)，在山峰地形的下底邊界位置(即地形突變的位置)處都會(huì)產(chǎn)生反映x、y方向邊界信息的正負(fù)極值虛假異常區(qū)域，Tzx的實(shí)、虛部在x方向上有對稱異常，實(shí)部和虛部正負(fù)異常出現(xiàn)相反的特征，這與山峰地形在橫向上空氣與地下介質(zhì)的邊界特征是相匹配的，與此同時(shí)，Tzy在y方向出現(xiàn)類似特征。因此，在進(jìn)行反演數(shù)值模擬時(shí)，需要考慮到地形因素產(chǎn)生的虛假異常對反演結(jié)果所產(chǎn)生的影響。

合成山谷純地形的算例進(jìn)行對比分析，其上頂面邊長為2.5 km，下底邊長為0.5 km，下底面低于地面0.45 km，其他計(jì)算條件及參數(shù)信息與山峰地形一致。如圖7所示，在地形突變的位置，出現(xiàn)了與山峰地形相反的ZTEM傾子平面響應(yīng)特征，此時(shí)Tzx和Tzy的實(shí)虛部的虛假異常區(qū)域要比山峰地形的更大，異常的極值也略大于山峰地形的虛假異常，這些虛假異常在起伏地形下的反演工作不可忽略。

圖6 頻率在50 Hz時(shí)山峰模型的ZTEM響應(yīng)Tzx和Tzy平面分布(白色框?yàn)槟Ｐ偷妮喞?Fig.6 Contour of the ZTEM tipper response Tzx & Tzy at 50 Hz generated from 3D in peak terrain model(the white box is the outline of the model)

圖7 頻率在50 Hz時(shí)山谷模型的ZTEM響應(yīng)Tzx和Tzy平面分布(白色框?yàn)槟Ｐ偷妮喞?Fig.7 Contour of the ZTEM tipper response Tzx & Tzy at 50 Hz generated from 3D in valley terrain model(the white box is the outline of the model)

3.2 反演算例分析

為了探討地形因素對反演結(jié)果可能造成的影響，分別設(shè)計(jì)一個(gè)山峰和山谷地形條件下的低阻單棱柱體模型進(jìn)行試算，從而進(jìn)一步說明所開發(fā)ZTEM反演代碼對帶地形條件下的地電模型反演的有效性和穩(wěn)定性。

模型一：圖8所示的是一個(gè)山峰地形條件下低阻棱柱體的反演算例，在山峰下方離水平參考地面0.2 km處有一個(gè)x×y×z為1 km×1 km×0.3 km的單棱柱體，其電阻率為100 Ω·m，地下介質(zhì)的電阻率為500 Ω·m，空氣層電阻率為108Ω·m。將研究區(qū)域離散化，目標(biāo)體中心區(qū)域x、y、z方向分別用250 m×250 m×100 m的網(wǎng)格進(jìn)行均勻剖分，由于遠(yuǎn)參考基站位置距離目標(biāo)體研究區(qū)域相對較遠(yuǎn)，因此，采取向外加大間距的方式進(jìn)行網(wǎng)格剖分，x、y方向均設(shè)置為28個(gè)網(wǎng)格，z方向有30個(gè)網(wǎng)格(包含11層空氣層)。

反演的初始模型為帶地形的均勻介質(zhì)空間，電阻率取500 Ω·m，使用頻點(diǎn)個(gè)數(shù)有5個(gè)(25、100、200、400、500 Hz)，水平磁場取值位置在(5 100 m，5 100 m，0 m)的水平地表處，垂直磁場取在離山峰頂面高100 m所處平面上，對所示模型進(jìn)行三維正演計(jì)算，模擬測點(diǎn)數(shù)有36個(gè)(圖8中倒三角形所示位置)，對參與反演的傾子分量加入5%的隨機(jī)誤差，用于模擬實(shí)測資料。經(jīng)過10次迭代后，歷時(shí)9 h40 min， 反演結(jié)果如圖9第一、第二行所示。

圖8 低阻棱柱體模型示意(山峰地形)Fig.8 Schematic diagram of conductive prism model (peak terrain)

第一行—5個(gè)頻率的水平地形ZTEM傾子反演結(jié)果；第二行—5個(gè)頻率的帶地形ZTEM傾子反演結(jié)果；第三行—8個(gè)頻率的帶地形ZTEM傾子反演結(jié)果；第一列—深為500 m的水平切片；第二列—x=0時(shí)沿y方向的垂直切片；第三列—y=0時(shí)沿x方向的垂直切片；黑色虛線—棱柱體的邊界the top row—the results from the inversion of tipper data (ZTEM) with five frequencies；the second row—the results from the inversion of tipper data (ZTEM) including topography with five frequencies；the third row—the results from the inversion of tipper data (ZTEM) including topography with eight frequencies；the first column—the horizontal slices at 500 m depth；the second column—the vertical slices at y=0 m along the y axis；the third column—the vertical slices at y=0 m along the x axis；the black dashed lines—the prism margins圖9 低阻體模型ZTEM傾子響應(yīng)數(shù)據(jù)不帶地形與帶地形的三維反演結(jié)果(山峰地形條件)Fig.9 Results from the 3D inversion of tipper ZTEM data generated from 3D conductive prism model (peak terrain)

從反演結(jié)果圖中可以發(fā)現(xiàn)，帶地形的ZTEM三維反演基本恢復(fù)出了目標(biāo)體的位置和形態(tài)，所恢復(fù)的模型在橫向邊界上具有較好的約束，縱向上存在拉伸的趨勢，對縱向邊界的恢復(fù)效果較差，說明ZTEM方法對縱向上的信息缺乏約束力。在山峰地形的影響下，由于在地形突變位置處存在虛假異常，因此實(shí)際恢復(fù)出的目標(biāo)體電阻率偏大(最小為114.7 Ω·m)。對于不帶地形的ZTEM三維反演，雖然基本恢復(fù)出了目標(biāo)體的主要特征，但是其位置和形態(tài)發(fā)生了較大畸變，反演效果較差，且縱向上的拉伸趨勢更為明顯，目標(biāo)體下方出現(xiàn)虛假異常，同時(shí)在地形起伏變化的位置附近產(chǎn)生了較為嚴(yán)重的虛假異常。

當(dāng)把頻點(diǎn)數(shù)增加至8個(gè)參與帶地形的反演而其他反演條件不變時(shí)，反演結(jié)果如圖9第三行所示，縱向上仍存在一定的拉伸，實(shí)際反演出的電阻率更接近目標(biāo)體電阻率，可達(dá)到95.2 Ω·m，對反演結(jié)果有一定的改善。

模型二：山谷地形條件下，距離地表0.65 km處有一個(gè)1 km×1 km×0.3 km的低阻單棱柱體，山谷底邊界距離地表有0.45 km，上頂邊長有2.5 km，如圖10所示，其他計(jì)算條件與山峰地形條件下的低阻單棱柱體一致。

經(jīng)過10次迭代，歷時(shí)12 h33 min，ZTEM傾子數(shù)據(jù)帶地形反演結(jié)果與不帶地形的反演結(jié)果對比如圖11所示。帶地形ZTEM三維反演所恢復(fù)的反演模型較好反映出了目標(biāo)體的形態(tài)與位置，尤其是對目標(biāo)體的橫向邊界位置恢復(fù)效果更為理想，相比于山峰地形條件，其縱向上的拉伸趨勢有所減弱，但對于上下底邊界的約束同樣存在不足，這是ZTEM方法本身所存在的一個(gè)缺陷。受山谷地形條件的影響，目標(biāo)體所恢復(fù)出來的電阻率離真實(shí)電阻率還有一定差距(最小為142.2 Ω·m)，這與山谷地形類似為高阻體的特征相符。

對于不帶地形的ZTEM三維反演模型,大致反演出了目標(biāo)體的輪廓，但目標(biāo)體形態(tài)和位置與正演模型相差較大，產(chǎn)生了較為嚴(yán)重的畸變，在目標(biāo)體周圍反演出了虛假異常體，并且在地形突變位置附近反演產(chǎn)生了部分虛假異常。

圖10 低阻棱柱體模型示意(山谷地形條件)Fig.10 Schematic diagram of conductive prism model (valley terrain)

第一行—水平地形ZTEM傾子反演結(jié)果；第二行—帶地形ZTEM傾子反演結(jié)果；第一列—深為650 m的水平切片；第二列—x=0時(shí)沿y方向的垂直切片；第三列—y=0時(shí)沿x方向的垂直切片；黑色虛線—棱柱體的邊界the top row—the results from the inversion of tipper data (ZTEM)；the second row—shows the results from the inversion of tipper data (ZTEM) including topography；the first column—the horizontal slices at 650 m depth；the second column—the vertical slices at x=0 m along the y axis；the third column—the vertical slices at y=0 m along the x axis；the black dashed lines—the prism margins圖11 低阻體模型ZTEM傾子響應(yīng)數(shù)據(jù)不帶地形與帶地形的三維反演結(jié)果(山谷地形條件)Fig.11 Results from the 3D inversion of tipper(ZTEM) data generated from 3D conductive prism model (valley terrain)

4 結(jié)論

基于三維有限差分正演算法分析了純起伏地形的ZTEM響應(yīng)特征，在考慮地形因素的影響下，開發(fā)實(shí)現(xiàn)了帶地形的三維ZTEM反演算法。合成算例分析表明，基于有限差分正演的帶地形ZTEM數(shù)據(jù)空間OCCAM反演方法可以得到比較接近地下真實(shí)導(dǎo)電性結(jié)構(gòu)的反演模型，說明了ZTEM帶地形進(jìn)行反演的算法準(zhǔn)確有效，可以有效減少虛假異常的產(chǎn)生。

通過山峰和山谷地形下低阻棱柱體的帶地形與不帶地形的反演算例計(jì)算結(jié)果對比分析，表明地形因素會(huì)對模型的形態(tài)、位置和電阻率恢復(fù)值帶來一定影響。水平地形的反演算法計(jì)算帶地形的數(shù)據(jù)時(shí)會(huì)帶來一些虛假異常，并使目標(biāo)體的形態(tài)產(chǎn)生畸變，而帶地形的三維ZTEM反演能夠得到比較接近真實(shí)的反演目標(biāo)體模型，可以有效減少虛假異常，尤其對模型橫向上的邊界具有較好的約束能力，但對垂向邊界的恢復(fù)能力略顯不足。此外，由于模擬起伏地形的需要，有限差分法需要在空氣與地表界面處做精細(xì)剖分，不可避免地會(huì)增加模型網(wǎng)格數(shù)，從而導(dǎo)致出現(xiàn)反演時(shí)間增長等問題。后續(xù)研究需要進(jìn)一步考慮采用并行算法來解決運(yùn)行效率的問題。

致謝感謝審稿專家提出的寶貴修改意見和編輯部的大力支持！