帶地形的ZTEM傾子資料三維正反演研究

李志強，孫洋，譚捍東，張承客

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

0 引言

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

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

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

1 起伏地形下ZTEM三維正演

1.1 傾子資料整理

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

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

(1)

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

(2)

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

(3)

1.2 傾子響應的計算

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

(4)

(5)

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

KE=S,

(6)

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

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法求解這兩個正演方程得到。

(8)

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

(9a)

(9b)

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

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

(9c)

1.3 正演精度的驗證

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

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

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

2 三維反演算法

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

(10)

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

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

(11)

(12)

圖3 二維山峰純地形模型傾子響應的二維、三維計算結果對比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反演問題中，為實現(xiàn)模型的更新迭代，需要通過計算得到的雅克比矩陣Jk獲得每次模型更新所需的模型修正量，而最為關鍵的是計算正演響應F和求解對應的Jk。傾子用于反演的第k個模型參數(shù)的偏導數(shù)Jk=(?F/?m)k可相應表示為Jk=(?T/?m)k[10]，將式(3)對模型電導率參數(shù)σk求偏導數(shù)，有:

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

(13)

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

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

(14)

(15)

接著將關系式(15)代入式(13)，則第k個模型電導率參數(shù)σk的擾動對第j個測點的改變量?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計算出后再求解

K·v1=1Gjzx,

(19)

K·v2=2Gjzx,

(20)

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

3 合成算例分析

3.1 三維純地形的異常特征

設計一個正四棱臺狀的山峰地形來分析其ZTEM正演傾子平面響應特征。如圖5所示，棱臺狀山峰的頂面相對于地面有0.45 km，上頂面邊長有0.5 km，下底邊長有2.5 km，空氣層電阻率為108Ω·m， 地下介質的電阻率為100 Ω·m， 選取地表固定點(5 750 m，5 750 m，0 m)作為兩個水平磁場分量取值位置，垂直磁場的取值高度設置在離山峰高100 m處。x、y、z方向網(wǎng)格單元剖分數(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三維傾子響應，圖中白色實線表示山峰下底邊界。由圖不難發(fā)現(xiàn)，在山峰地形的下底邊界位置(即地形突變的位置)處都會產生反映x、y方向邊界信息的正負極值虛假異常區(qū)域，Tzx的實、虛部在x方向上有對稱異常，實部和虛部正負異常出現(xiàn)相反的特征，這與山峰地形在橫向上空氣與地下介質的邊界特征是相匹配的，與此同時，Tzy在y方向出現(xiàn)類似特征。因此，在進行反演數(shù)值模擬時，需要考慮到地形因素產生的虛假異常對反演結果所產生的影響。

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

圖6 頻率在50 Hz時山峰模型的ZTEM響應Tzx和Tzy平面分布(白色框為模型的輪廓)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時山谷模型的ZTEM響應Tzx和Tzy平面分布(白色框為模型的輪廓)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 反演算例分析

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

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

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

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

第一行—5個頻率的水平地形ZTEM傾子反演結果；第二行—5個頻率的帶地形ZTEM傾子反演結果；第三行—8個頻率的帶地形ZTEM傾子反演結果；第一列—深為500 m的水平切片；第二列—x=0時沿y方向的垂直切片；第三列—y=0時沿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傾子響應數(shù)據(jù)不帶地形與帶地形的三維反演結果(山峰地形條件)Fig.9 Results from the 3D inversion of tipper ZTEM data generated from 3D conductive prism model (peak terrain)

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

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

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

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

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

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

第一行—水平地形ZTEM傾子反演結果；第二行—帶地形ZTEM傾子反演結果；第一列—深為650 m的水平切片；第二列—x=0時沿y方向的垂直切片；第三列—y=0時沿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傾子響應數(shù)據(jù)不帶地形與帶地形的三維反演結果(山谷地形條件)Fig.11 Results from the 3D inversion of tipper(ZTEM) data generated from 3D conductive prism model (valley terrain)

4 結論

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

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

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